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Summary of Changes to this Version 

The last version of this document (4.2, October 2002) was written prior to the launch of ICESat and 
the analysis algorithms described therein were based on theory and tested with simulated data. 
After launch and the acquisition of real data, many changes and additions were made to the 
algorithms, both in terms of increasing the accuracy of computed parameters and the addition of 
new parameters. With regard to the latter, the laser problems encountered during the mission (loss 
of 532 nm laser energy) required the addition of codes to obtain as much information as possible 
from the 1064 channel. Originally, all GLAS atmospheric data products were to be obtained from 
the 532 nm channel only. Early in the mission it became apparent that the 1 064 data would have to 
be used if any substantial quantity of atmospheric data were to be obtained by the mission. Many 
of the additions to this version are related to the use of the 1064 channel to retrieve atmospheric 
quantities. However, the signal quality of the 1064 channel limits the products that can be obtained 
to cloud height, relatively thick aerosol layer heights, blowing snow detection and total column 
optical depth. Boundary layer height and optical properties of clouds and aerosols require the 
higher signal to noise ratio afforded by a nominally functioning 532 channel, which unfortunately, 
existed for only a short time. 

While there were numerous changes made to the codes since October of 2002, the major changes 
to this version can be grouped into the following areas: 

1 ) Addition of 1 064 channel derived products (Section 3.3.2) 

2) Cloud /Aerosol discrimination (Section 3.3.1) 

3) 532 channel background computation (Section 3.1.1) 

4) 1064 channel droop correction (Section 3.1.1) 

5) Extinction calculation (Section 3.6) 

6) Blowing Snow detection (Section 3.5) 
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1 1ntroduction 


Launched in January 2003, the Geoscience Laser Altimeter System (GLAS) is the only instrument 
aboard the ICESat satellite and is an atmospheric lidar in addition to a surface altimeter. GLAS 
operated roughly 3 times per year in month-long periods from February, 2003 to October, 2009 
providing high resolution measurements of global topography with special emphasis on the 
determination of the temporal changes of ice sheet mass over Antarctica and Greenland. The 
primary atmospheric science goal of the GLAS cloud and aerosol measurement is to determine the 
radiative forcing and vertically resolved atmospheric heating rate due to cloud and aerosol by 
directly observing the vertical structure and magnitude of cloud and aerosol parameters that are 
important for the radiative balance of the earth-atmosphere system, but which are ambiguous or 
impossible to obtain from existing or planned passive remote sensors. A further goal is to directly 
measure the height of atmospheric transition layers (inversions) which are important for dynamics 
and mixing, the planetary boundary layer and lifting condensation level. Towards these goals, the 
various level 1 and 2 atmospheric data products which are generated on the Investigator-led 
Science Information Processing System (ISIPS) are: 

1 . GLA02 - Normalized relative backscatter (1064 and 532) 

2. GLA07 - Calibrated attenuated backscatter cross section (1064 and 532) 

3. GLA08 - Planetary Boundary Layer (PBL) height and elevated tropospheric aerosol layer height 
(derived from the 532 channel) 

4. GLA09 - Cloud top (and bottom when possible) heights (1064 and 532) 

5. GLA10 - Attenuation-corrected backscatter and extinction cross section (532 only) 

6. GLA11 - Thin cloud and aerosol layer optical depth and total column optical depth (532 and 
1064) 

Because of the laser issues, GLAS was not operated continuously, but rather obtained data 3 
times per year in 33 day long observation periods. During these six years and all 16 operational 
periods (LI A, L2A - L2D, L3A - L3K), GLAS obtained high quality altimetry measurements with 
only minor loss of data as laser energy decreased. The atmospheric measurement, however, 
requires more laser energy and the quality of these measurements decreased considerably with 
time. This was especially true of the 532 channel which operated at full or near-full signal strength 
only for laser operation period L2A (October - November, 2003) and the first 2 weeks of the L2B 
campaign. The 532 channel is used to obtain cloud height, boundary layer height, attenuation 
corrected backscatter coefficient, aerosol and cloud layer optical depth, and extinction profiles. 
After the L2B observation period (17 Feb 04 - 21 Mar 04) the quality of the 532 nm signal is such 
that the 532 nm derived data products can only be gen erated from nighttime data. After 
observation period L3E, no 532 based data products are available. The 1064 laser energy did not 
decrease as rapidly and the data products that are derived from this channel maintained better 
consistency. However, the 1064 nm atmospheric products are limited to cloud layer height, aerosol 
layer height and 1064 total column optical depth (over oceans and ice sheets only). Extinction 
profiles are not generated from the 1064 channel. These products maintain a reasonable quality up 
through observation period L3I (02 Oct 07 - 05 Nov 07). Table 1 lists the dates of the 16 
observation periods and gives a qualitative assessment of the data quality. 
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After a short introduction, we will provide an overview of the GLAS instrument and prior lidar work 
which is pertinent to the GLAS data products discussed here, before presenting the details of the 
individual algorithms in section 3. Section 4 di scusses the practical applications and 
implementation issues of each algorithm including examples of output. Section 5 discusses ranging 
errors caused by multiple scattering of laser photons as they travel through the atmosphere. 


Table 1. GLAS Observa 

tion Periods and atmospheric data quality 

Date 

532 nm Channel 

1064 nm Channel 

Obs Period 

20Feb03 - 29Mar03 

None 

Excellent 

1 

25Sep03 - 18Nov03 

Excellent 

Excellent 

2A 

17Feb04 - 21Mar04 

Excellent - Good - Fair 

Excellent - Good 

2B 

18May04 - 21 Jun04 

Fair - Night only 

Marginal 

2C 

040ct04 - 09Nov04 

Fair - Night only 

Excellent 

3A 

17Feb05 - 24Mar05 

Fair - Night only 

Excellent 

3B 

20May05 - 24Jun05 

Fair - Night only 

Excellent 

3C 

21Oct05 - 24Nov05 

Fair - Night only 

Good 

3D 

22Feb06 - 28Mar06 

Fair - Night only 

Good 

3E 

24May06 - 26Jun06 

None 

Fair 

3F 

25Oct06 - 27Nov06 

None 

Fair 

3G 

12Mar07 - 14Apr07 

None 

Fair 

3H 

020ct07 - 05Nov07 

None 

Fair 

31 

17Feb08-21Mar08 

None 

Marginal 

3J 

040ct08 - 19Oct08 

None 

Poor 

3K 

24Nov08 - 17Dec08 

None 

Poor 

2D 








2 Overview and Background 

2.1 History 


The purpose of this document is to present a detailed description of the algorithm theoretical basis 
for each of the GLAS data products. This version (V5.0) is a I ong overdue update to the last 
version (V4.2) of the Atmospheric product ATBD which was completed prior to launch in 2003. This 
will be the final version of this document. The algorithms were initially designed and written based 
on the authors’ prior experience with high altitude lidar data on systems such as the Cloud and 
Aerosol Lidar System (CALS) and the Cloud Physics Lidar (CPL), both of which fly on the NASA 
ER-2 high altitude aircraft. These lidar systems have been employed in many field experiments 
around the world and algorithms have been developed to analyze these data for a n umber of 
atmospheric parameters. CALS data have been analyzed for cloud top height, thin cloud optical 
depth, cirrus cloud emittance (Spinhirne and Hart, 1990) and boundary layer depth (Palm and 
Spinhirne, 1987, 1998). The successor to CALS, the CPL, has also been extensively deployed in 
field missions since 2000 including the validation of GLAS and CALIPSO. The CALS and early 
CPL data sets also served as the basis for the construction of simulated GLAS data sets which 
were then used to develop and test the GLAS analysis algorithms. 

After launch in 2003, there were numerous updates, additions and fixes to the atmospheric data 
product codes which were then based on the GLAS data itself. Many of these changes were minor, 
such as the selection of better threshold values for layer detection and cloud aerosol 
discrimination. However, some were quite major like what was done to correct for a range 
dependent background in the 532 channel (see section 3.1. 1.2), blowing snow detection (section 
3.5) and the addition of 1064 derived products. After each software update and prior to its release, 
the codes were extensively tested and the output compared with observations (when available) 
and the prior version to validate the effectiveness of the changes while maintaining consistency. In 
all there were 33 versions of the software released with version 12 being the first version released 
after launch. 

2.2 Instrument and Data Description 

The GLAS atmospheric measurements are obtained from the 600 km polar orbiting platform both 
day and night using two separate channels. The 532 nm, photon counting channel is the most 
sensitive and provides the highest quality data when the laser energy is above about 10 mJ. 
Unfortunately, this level of laser energy was maintained only during the L2A and first half of the 
L2B observational campaigns. This channel employs an etalon filter which is actively tuned to the 
laser frequency, providing a very tight bandpass filter of about 30 picometers. This, together with a 
very narrow (180 fxr) receiver field of view (FOV), enables high quality daytime measurements 
even over bright background scenes. There are 8 separate photon counting detectors on GLAS, 
but 4 of these detectors failed during ground vibration testing. The return signal is split equally 
between these detectors. Thus, in addition to the lower than anticipated 532 nm laser energy, half 
of the 532 nm return signal is essentially discarded. Even with these problems, the 532 channel 
provided good data (capable of the retrieval of all atmospheric parameters) both day and night 
down to about 10 mJ laser energy. Below this point, the atmospheric retrievals cannot be reliably 
performed for daytime data. The nighttime retrievals are good down to about 4-5 mJ laser energy. 
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The 1064 nm channel uses an Avalanche Photo Diode (APD) detector with a much wider (0.1 nm) 
bandpass filter and FOV (475 pir). The sensitivity of the 1064 channel is limited by the inherent 
detector noise. However, experience has shown that the 1064 channel is able to measure 
backscatter cross section down to about 2.0x1 O' 6 nr 1 sr 1 , which means that it can detect fairly thin 
clouds and moderately thick aerosols down to an optical depth of about 0.10 for a 1.5 km thick 
layer. Table 2 lists the major GLAS system parameters which ultimately affect system performance 
and data quality. Note that the laser energy shown in Table 2 refers to the start of the mission. 
Significant laser energy reduction occurred with time for each of the three lasers. Also note that for 
laser one, which failed after 5 weeks of operation in March, 2003, the 532 channel detectors were 
not powered on because it was feared that outgassing from adhesives could potentially cause 
harm to the detectors. 

The GLAS laser transmits short (5 nanosecond) pulses of laser light (in the nadir direction) 
producing a footprint 70 meters wide upon striking the surface, and each footprint is about 175 
meters apart. The backscattered light from atmospheric clouds, aerosols and molecules is digitized 
at 1.953 MHz, yielding a vertical resolution of 76.8 meters. The horizontal resolution is a function of 
height. For the lowest 10 km, each backscattered laser pulse is stored. Between 10 and 20 km, 8 
shots are summed, producing a horizontal resolution of 5Hz or 1.4 kilometers. For the upper half of 
the profile (20-40 km), which is entirely within the stratosphere, 40 shots are summed, providing a 
horizontal resolution of about 7.5 kilometers. This approach was adopted for a number of reasons. 
First, the atmospheric processes of interest have more variability and smaller scales in the lower 
troposphere (particularly the boundary layer) than in the mid and upper troposphere. Second, the 
amount of molecular and aerosol scattering in the upper troposphere and stratosphere is so small 
that summing multiple shots is required to obtain a non-zero result. Lastly, this approach helps to 
reduce the amount of data that has to be stored on board the spacecraft and transmitted to the 
ground. 


Table 2. GLAS System Parameters 


Parameter 

532 Channel 

1064 Channel 

Orbit Altitude 

600 km 

600 km 

Laser Energy 

25 mJ 

70 mJ 

Laser Divergence 

1 1 0 jli rad 

1 1 0 |ii rad 

Laser Repetition Rate 

40 Hz 

40 Hz 

Effective Telescope Diameter 

95 cm 

95 cm 

Receiver Field of View 

1 80 |ii rad 

475 |ii rad 

Detector Quantum Efficiency 

60% 

35% 

Detector Dark Current 

3.0x1 0 -16 A 

50.0x1 0- 12 A 

RMS Detector Noise 

0.0 

2.0x1 0- 11 

Electrical Bandwidth 

1. 953x1 0 6 

1.953x106 

Optical Filter Bandwidth 

0.030 nm 

0.800 nm 

Total Optical Transmission 

30% 

55% 
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2.3 Description of GLAS Atmospheric Channel Data 

The atmospheric channel of GLAS provides a record of the vertical structure of backscatter 
intensity from the ground to a height of about 41 km with 76.8 meter vertical resolution. Two 
channels are employed, the Nd:Yag fundamental wavelength of 1064 nm and the frequency 
doubled 532 nm wavelength in the visible portion of the spectrum (green channel). The basic 
equation which describes the atmospheric return signal p(z) is the standard lidar equation 


(2.1) p(z) 


CEp(z)T 2 (z) 

2 ^ Pb ^ Pd 


where (3(z) is the total atmospheric backscatter cross section at an altitude z, T(z) is the 
transmission from the top of the atmosphere to altitude z, r is the range from the spacecraft to the 
altitude z, E is the transmitted laser pulse energy and C is a dimensional constant referred to as 
the lidar calibration constant. There are two range independent background terms, pb from 
scattered solar radiation and pd for any detector dark signal or noise. In the case where p would be 
the signal in watts returned to the receiver detector, the calibration constant is given as 

(2.2) C=cAT s /2 

where c is the light speed constant , A the area of the receiver and T s the optical transmission of 
the receiver system. For the GLAS 532 nm atmospheric channel the signal will be acquired as the 
photo-electron count rate from the detector n(z). In this case the calibration constant is given as 

(2.3) C=AT s lq/2h 

where X is the wavelength, q is the photon detection probability or quantum efficiency, and h is the 
Plank constant. The background radiance signal in terms of photo-electron count rate will be 

(2.4) nb=AT s lbQA/hc 

where lb is the background radiance and Q is the receiver solid angle and A is the optical 
bandwidth. The additional background signal will be any detector dark photo-electron count rate 
nd. 

The 1064 urn detector for GLAS is the same silicon APD detector that will be used for the surface 
return signal although a separate lower speed A/D signal acquisition will be used. The signal in 
this case is a voltage from the detector amplifier V(z). The calibration constant will be 

(2.5) C=AT s crg v /2 

where r is the detector responsivity in amps/watt, and g v is the voltage gain of the detector 
preamplifier. The detector background signal will be idg v where id is the detector dark current. 

The accuracy of the received GLAS atmospheric signals is limited by the fundamental probability, 
or signal shot noise of the signal. For the case of the 532 nm photon counting signal, the noise 
factor is given by Poisson statistic. The signal to noise ratio will then be given by 
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(2.6) S / N = 


n(z) 


V n ( z ) + n b +n d 


Where n(z) is the number of photons detected by the lidar at range z, nb is the background signal 
and nd id the detector dark count (noise). In the case where the signal is voltage derived from a 
detected current the basic signal to noise will be: 


S/N= , 

(2J) V 2Af ( 1 s+ib+i d )e 

where i s is the detector current produced by the backscattered signal, ib is the detector current 
produced by background ambient light id is the detector dark current, A/ is the system electronic 
bandwidth, and e is electron charge. The signal noise defines the degree to which the lidar data 
may be usefully applied. 



3 GLAS Atmospheric Algorithms 

This section will address in detail the structure and content of the six algorithms which 
comprise the level 1 A, 1 B and level 2 GLAS atmospheric data products. A theoretical 
description will be given for each algorithm followed by error quantification and a 
description of the confidence (quality) flags which attempt to assign a confidence level 
to the quality of the algorithm output. Section 4 will discuss the issues related to the 
practical application and implementation of the algorithms. 

3.1 Normalized Lidar Signal (GLA02) 

3. 1. 1 Theoretical Description 

3. 1. 1. 1 Normalized Lidar Signal - LI A 

The normalized lidar signal is a level 1A data product which applies the fundamental 
corrections and normalizations to the raw data as well as providing an estimate of the 
height of the first cloud top and/or the bin location of the ground return. Additionally, it 
flags each 532 nm channel bin which has reached saturation so that it may be 
corrected in later processing. The algorithm applies range and laser energy 
normalizations, computes and subtracts out the ambient background signal, and 
performs dead time correction to the photon counting (532 nm) channel. The dead time 
correction is performed by using 8 separate and unique look-up tables which contain a 
dead time corrected value for each possible output from the photon counting detector. 
The dead time look-up table was constructed using manufacturers laboratory 
measurements of the performance of each detector. In the case of the 1064 channel, 
the digital counts that are output from the analog to digital converter must first be 
converted back to a voltage using a lookup table which has been calibrated and tested 
in the laboratory. The background subtraction, energy and range corrections are then 
applied to the data. 
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The basic output of GLA02 is the generation of what we call normalized lidar signal (P'(z)). From 

( 2 . 1 ) we first subtract the background, then multiply by the square of the range (in km) from the 
lidar receiver to the return bin (R 2 ) and divide by the laser energy (E, in millijoules). Here, we have 
combined the detector dark current (Pd) and the ambient background light (Pb) into one background 
term (B). We must also perform dead time correction on the raw photon counts (for the 532 
channel) and convert from digital counts to volts in the case of the 1064 channel. Now, a further 
consideration for the 532 channel is the etalon transmission. For the 532 channel, a narrow-band 
etalon filter is used for rejection of background light. The etalon bandpass is about 30 picometers 
wide. The laser frequency may shift considerably on a shot to shot basis, which could result in a 
loss of return signal, since the laser frequency is not centered on the fliter bandpass. On board the 
spacecraft, a part of the laser energy will be split off and sent through the etalon. The amount of 
energy passing through the etalon will be m easured and s ent down in the telemetry. In the 
telemetry spreadsheet this is known as “Dual Pin A.” The ratio of this to the outgoing laser energy 
at 532 times a calibration coefficient gives us a relative measure of the etalon transmission. The 
calibration coefficient (y) will be determined by laboratory measurement. Thus, if we let a = y DPA / 
E532 where DPA is the “Dual Pin A” output, the equation to produce the normalized signal for the 
532 channel is: 

(3.1.1) p; 32 (z) = C 532 /J 532 (z)T 2 s 32 = ( DC[S 532 (Z )] - DC[B(z) 532 ])R 2 /( aE 532 ) 

where DC denotes the dead time correction lookup table discussed above. Note that in the 
denominator, the transmit energy cancels as: aEs 32 =yDPA 

For the 1064 channel, we must first convert the digital counts to voltage for both the background 
(Vb) and the atmospheric return signal (V s ) before computing the normalized signal. 

3.1.2 V b =(B m4 F v -V 0 )/(GA) 

3.1.3 V s (z) 

— (^1064 ( Z )F V -V 0 )/(GA) 

where B 1064 is the 1064 nm background (computed from equation 3.1.6 below), F v is a constant 
(0.01560) relating digital counts to volts (volts per count), Vo is the voltage offset (currently set to 
0.90), G is the amplifier gain (currently set to 18.0) and A is the 1064 programmable attenuation 
setting (which have values of: 1, 1/1.77, 1/3.16, 1/5.6, and 1/10). Next we compute the normalized 
return signal 

(3.1 .4) P; 064 (z) = C m4 j3 m4 (z)T 2 io 64 = D r (V s (z) - V b )R 2 / E l064 


The range from the spacecraft to the return bin (R) should be in kilometers and the laser energy (E) 
is in millijoules. The voltage must then be multiplied by the detector responsivity factor (D r = 
4.4x1 0 7 ) which has units watts per volt. The units for the 1064 channel are watts*km 2 /mJ. The units 
for the 532 channel are photons/bin*km 2 /mJ. 
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Note also that the 8 shot and 40 shot summed data from the middle and upper layers must be 
normalized by the number of shots summed in that layer before equations 3.1.1 and 3.1.4 can be 
applied. Thus, all data from the 5 Hz middle layer must be divided by 8 and all data from the upper 
layer (532 only) must be divided by 40 before application of equations 3.1 .1 and 3.1 .4. 

3. 1.1. 2 Background Computation 

The background signal (B) for the two channels is computed for each laser shot (40 Hz) from time 
integrated measurements of the background intensity at two separate times relative to laser fire. 
The first is prior to the laser beam reaching the atmosphere (about 100 km altitude), and the 
second is well after the beam strikes the Earth (approximately 40 km below ground). The two 
background measurements for each channel are stored as two byte values and must be 
normalized before use in equations 3.1.1 and 3.1.4. Letting Tb equal the background integration 
time in microseconds, and I 532 and I 1064 the integrated background signal for Tb microseconds, the 
background values (counts per bin) to be used in equations 3.1 .1 and 3.1 .4 are: 

(3.1.5) B 532 =I 532 /(l.953T b (A)) 

-6) 064 = A 064 /(l.953T b (A)) 

The background integration time is not the same for the two channels. Tb(532) is 256 
microseconds, while Tb(1064) is half that or 128 microseconds. The background is computed in this 
manner for the two integration periods. Originally, it was thought that an average of the upper and 
lower backgrounds would be used for the 532 channel, but after launch it was discovered that the 
gating off and then back on of the 532 SPCM detectors (between shots) induced a time dependent 
responsivity to the detectors that was caused by heating of the detectors from the incident 
background photons (this problem only exists for the 532 channel). The detector responsivity 
change is a function of the background level itself. It causes for instance, the upper background to 
be greater than the lower background for a certain range of background level but for another range 
of background level, the upper background would be lower in value than the lower background 
value. This made it impossible to utilize the background measurements either individually or as an 
average of the two. A scheme was adopted to address this by trying to compute the time (range) 
dependent responsivity of the 532 channel detectors. In essence, we compute a range dependent 
background profile comprised of bins of the same physical dimension as the signal bins of the 
atmospheric backscatter profile. Each bin of the signal profile has a corresponding background bin 
which is subtracted from the signal. Two schemes are used to compute the time (range) dependent 
background. One is a linear fit between the upper and lower backgrounds and the second utilizes a 
third background point computed from the upper 10 km of the 41 km 532 channel signal profile (31 
to 41 km altitude). A polynomial is then fit to the three points to obtain the range dependent 
background profile. Note that this background profile is not stored on the product. 

Though much effort was put into trying to completely solve this background problem, for certain 
lighting conditions it became apparent that neither of the above approaches were able to entirely 
eliminate the problem. The most troublesome condition occurs when the spacecraft goes from 
looking at a relatively low albedo surface to a high albedo surface very quickly. Also, as the 
spacecraft emerges from darkness into daylight over the polar regions is another difficult situation 
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where this background problem tends to be most acute. The overall effect of not being able to 
completely compensate or correct for it is that the 532 signal may at times not be well calibrated 
over areas of highly variable background. 

The computed signal profiles defined by 3.1.1 and 3.1.4 have the same horizontal resolution as the 
raw input data. This means that from -1 to 10 km altitude, both the 532 and 1064 channels will be 
40 Hz, between 10 and 20 km the profiles will be at 5 Hz, and between 20 and 41 km we have only 
532 data at 1 Hz. Note that the background computation as described above will be performed at 
40, 5 and 1 Hz (the 5 and 1 Hz backgrounds are computed by averaging the 40 Hz background 
measurements) for the 532 channel and at 40 and 5 Hz for the 1064 channel. Care must be taken 
to use the appropriate background in equations 3.1.1 and 3.1.4 as dictated by the given layer (40 
Hz background for the lowest layer, 5 Hz background for the middle layer and 1 Hz background for 
the upper layer). This also applies to the laser energy as well, as it is reported at 40 Hz. The 5 and 
1 Hz laser energies must be calculated and the appropriate one applied according to which layer is 
being processed. Both the background and laser energies at 40, 5 and 1 Hz are stored as part of 
the GLA02 output. 

In addition to the background being calculated from the high and low integration periods, it is also 
calculated from the last 8 bins of the lidar profile for both 532 and 1 064. A fourth element is added 
to the background array stored on the product: 

BG(1)- upper background 

BG(2) - lower background 

BG(3) - Background used in computing NRB 

BG(4) - Background computed from the last 8 bins of the profile 

Note that when the flag is set to compute and use a range dependent background as described 
above, BG(3) as defined above does not contain the background used in computing the 532 
channel NRB and this range dependent background is not stored on the GLA02 product. The 1064 
channel does not have this range dependent background problem and the background elements 
stored on the products are as defined above. However, since the 1064 channel is AC coupled, the 
background is electronically removed and instead a constant offset (value = 54.47) is electronically 
added. Thus, the 1064 background as defined by equation 3.1.6 is not used. Rather the constant 
value of 54.47 is used for the background value in equation 3.1.4. The 1064 channel, however, has 
its own interesting problem that was discovered after launch. This is described in the next section. 

3.1. 1.3 1064 Channel Droop Correction 

The 1064 channel is AC coupled and suffered from an effect that became known as signal droop. 
In essence, sometime after the detector is hit by a substantial signal from clouds it will start to lose 
signal after a certain time so that the signal does not return to the zero level but instead goes 
considerably below that for a fairly large time (10’s of microseconds). Figure 3.1.1a shows an 
example of this effect. The large signal at about 7-8 km altitude is caused by a relatively thick 
cloud. Note how after the maximum signal is attained, it rapidly drops off to a point below zero for a 
vertical distance of about 6 km in this case. Finally, near 1 km altitude, the signal has recovered to 
a near zero value. The vertical distance (or time) it takes to recover is dependent on the initial 
signal strength that hits the detector. This in turn is related to the cloud optical depth. The droop 
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effect is due to an ill-designed electronics circuit in the 1064 detector chain. In the GLAS LI A 
(GLA02) processing a correction for this effect has been coded but it does not completely recover 
the lost signal. The degree to which the signal is recovered depends on the magnitude of the initial 
signal on the detector. For medium to thin optical depth clouds (like most cirrus), the correction 
works very well. For thick clouds like the one shown in figure 3.1.1, the signal is only partially 
corrected. This can be seen in Figure 3.1.1b which is the output of the correction routine for the 
raw signal shown in Figure 3.1.1a. 



Figure 3.1.1. Example of the 1064 nm channel response after the detector is hit by a large signal 
from a thick cloud showing signal droop (a) and the degree to which the droop correction algorithm 
is able to correct this effect (b). For a complete correction, the signal below cloud bottom (6.5 km 
height) in (b) would be very near zero, as it is above 9 km. 

The correction routine first finds the first cloud return in the profile. Call this bin kk. It then subtracts 
off the baseline for the 1064 channel which was determined from laboratory measurement to be 
54.47 and performs a double integration of the signal as shown in the code snippet below. Here, 
signal_1 064 is the raw digitized counts (0-255) of the 1064 channel. 

d_raw_signal = signaM 064 - 54.47 ; 54.47 is the zero signal level of the A/D output 

d_integ_signal = O.OdO 
do i_pass = 1 ,2 
d_sumd = O.OdO 
do k=kk,gi_20_m1km 
d_sumd = d_sumd + d_raw_signal(k) 
d_integ_signal(k) = d_sumd 
enddo 

d_integ_factor = 9.5d0 / 325.0d0 
if (i_pass == 2) d_integ_factor = d _integ_factor / 9.0d0 
d_integ_signal = d_integ_signal * d_integ_factor 
do k=kk,gi_20_m1km 

d_cor_signal(k) = d_raw_signal(k) + d_integ_signal(k) 
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enddo 

d_raw_signal = d_cor_signal 
enddo 

signal_1064 = d_raw_signal + 54.47 

The droop correction is performed before background subtraction, energy normalization and range 
correction (ie. before application of equation 3.1.4). 

3.1. 1.4 532 nm Channel Afterpulse 

Another problem that was discovered after launch affected the 532 channel. After hitting very 
dense clouds and the ground, sometimes the 532 detectors become “stuck on” and put out a high, 
continuous signal level for a few microseconds. This is known as detector afterpulsing. This 
problem only happened for a very large signal from a thick cloud or the ground and only affects 
data bins that would normally be totally attenuated. Thus, in essence, this problem is really only 
cosmetic as the bins affected would not contain useful signal. The corrective code, implemented in 
the LAI processing, identifies the affected bins and zeroes them out. Recognizing when this effect 
occurs is not difficult as the detectors return to normal operation (signal output) not all at once, but 
in a series of discrete steps each lasting a few microseconds. The effect does occur from thick 
clouds but is more often produced from the ground returns over highly sloping terrain. An example 
of the afterpulse effect on the raw data is shown in Figure 3.1 ,2a and the output from the correction 
procedure is shown in Figure 3.1.2b. 



Figure 3.1.2 An example of the 532 nm detector afterpulsing effect when encountering thick clouds 
and (usually only) sloping terrain, (a) shows the raw data and (b) the corrected data. 


3.1.1.5 Calibration Pre-Processing, Predicted Cloud Height and Ground Return Bin 

Accurate calibration of the lidar channels is a v ery important step that is required to obtain 
quantitative information on clouds and aerosols such as optical depth or extinction. The 532 
channel can be calibrated using the molecular backscatter signal at that wavelength in particulate- 
free regions of the atmosphere. The 1064 channel cannot use this method because the 1064 nm 
molecular backscatter cross section is too small to measure. The 1064 calibration constant is 
computed from validation campaigns utilizing simultaneous, co-located measurements from the 
Cloud Physics Lidar (CPL) on the NASA ER-2 aircraft. It is also computed using calibrated 
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measurements of cirrus clouds from the 532 channel with the assumption that the backscatter 
cross section of cirrus clouds is wavelength independent. For the 532 channel, in the calibration 
pre-processing step, we calculate the average normalized lidar signal at two calibration heights for 
segments of data roughly 10 minutes in length continuously throughout the entire orbit. The first 
calibration height is constant at around 25 km (read in from the constants file) and will be used for 
the 532 channel only. The second calibration height will be calculated from the data comprising 
that segment and will occur at the height of the minimum average signal between 8 and 15 km. 
The calibration pre-processing described below is not a part of the GLA02 process, but instead is 
implemented as a stand alone module that runs after GLA02 completes. It uses the output of 
GLA02, namely the normalized signal discussed in section 3.1 .1 For the purpose of discussion, we 
will call the calibration pre-processing module the ‘Segment Averaging Module’ or SAM for short. 
After SAM is run, another module is run which uses the output of SAM to compute the actual 
calibration constants. The general processing scenario for SAM is as follows: 

1 . Construct a 1 Hz continuous profile of P’ from -1 to 41 km for the 532 channel and from -1 to 
20 km for the 1 064 channel. 

2. Add the background to ‘summing’ variables for each channel 

3. Sum the P'532 data from Hi to H2 km and add it to a ‘summing’ variable. The values of Hi and 
H 2 will be roughly 24 and 27, respectively, but is changeable and read in from the constants 
file. Increment an ‘upper counter’. 

4. Check for clouds from 22 km to 8 km above ground. If clouds were not found for this second, 
then do the following (number 5 below): 

5. Add the 1 Hz data (each bin) between 8 and 15 km to a 'summing' array for each channel. 
Increment a ‘ lower counter’. 

6. If you have been doing this for t minutes, where t is read in from the constants file (default 
value: t=1 0), and at least 50 percent of the expected number of seconds have been summed 
(based on the 'upper counter’), then do the following: 

a. compute the average 532 signal from Hi to H 2 km for the entire ‘t’ minute segment. Call 
this P2(532) from the sum generated in step 3 above. 

b. If the ‘lower counter’ exceeds 50 percent of the expected number of seconds, then perform 
c,d, and e below. Otherwise, set PI (532) and PI (1064) to invalid and skip c,d and e. This 
effectively means that clouds have made calculations impossible at the lower height. 

c. Compute the average 532 and 1064 profiles between 15 and 8 km from the summing array 
produced in steps 4 and 5 above. 

d. Find the height of the minimum in the 532 average profile between 8 and 15 km call this 
hmin - this is the lower calibration height 

e. Compute the average of the data between hmin+D and hmin-D km for both the 532 and 
1064 channels, where D is in km and is read from the constants file (default = 1km). Call 
these P1(532) and P1(1064). 

f. Compute the average background for the segment for each channel call these B532 and 
B1064 

g. Output to a file: P1(532), P1(1064), P2(532), B532, B1064, hmin, D, Hi, H 2 and: the 
latitude, longitude and time at ‘m’ points along the segment, where m is a variable read 
from the constants file, not to exceed 30. A default value for m is 20. These points would 
be t/m minutes apart. 
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7. If after ‘t’ minutes, less than 50 percent of the expected number of seconds have been 
summed (based on the ‘upper counter’), then output missing values (invalid) for PI (532), 
PI (1064), P2(532), B532, B1064, and the other output described in 6g above. 

8 Zero out summing variables, summing array and counters 

9 Process next ‘t’ minute segment in the same manner 

(3.1 .7) PC = H sat - [H mm -H offmm ] -41 km + PC bias 

A complication arises in that the data read in by GLA02 are not vertically aligned from second to 
second. Onboard the spacecraft, the start of data (the height above the ellipsoid of the top most 
bin) is calculated from equation 3.1.7. The spacecraft position (updated every second) is used to 
retrieve the Digital Elevation Model (DEM) value corresponding to the spacecraft location. The 
DEM is a 1x1 degree land surface elevation in km above the ellipsoid. Hmin - Hoftmin represents the 
minimum elevation for a particular grid box, and H sa t is the height of the spacecraft above the 
ellipsoid from onboard GPS measurements. The subtraction of 41 km in equation 3.1.7 insures that 
the top of the profile for any given second will be 41 km above the minimum ground elevation for 
the current DEM grid box. This also means that the bottom bin of the profile will be 1 km below the 
minimum elevation for the current DEM grid box (since the total lidar profile is 42 km in length). 
PCbias can be used to shift the whole lidar profile up (when negative) or down (when positive), but 
will normally be zero. Because this process is happening each second, the bin that corresponds to 
a given height above mean sea level may change from second to second. Thus, to accomplish 
steps 4 and 5 above, the DEM value that was used onboard the spacecraft must be known to SAM 
so that it can compute the height of a given lidar bin number. Either that or simply the range from 
the spacecraft to the start of data. For the case where the onboard DEM value is used, H = (548- 
n)*0.0768 + [Hmin - Hoffmin] - 1.0 - PCbias, gives the rough height in km above the ellipsoid for any 
bin n, where n=1 is the topmost lidar bin of the complete profile. 548 represents the total number of 
bins in a complete 532 nm profile (as constructed in step 1 above). Thus, as an example, the bin 
number corresponding to H = 36 km would be: n = 548 - ((36 + 1.0 + PCbias + Hoffmin “ Hmin) / 
0.0768). Note that these heights are calculated with respect to the ellipsoid, which can depart as 
much as 200 meters from mean sea level. 

The intent of this process (SAM) is to produce the average signal (P’ 532 ) at the two calibration 
heights and PW at the lower calibration height every ‘t’ minutes. Depending on the magnitude of 
‘t’, this will correspond to about 6-10 points per orbit. The file created by SAM will be read in by a 
‘CALibration Module’ (CALM) that will produce the calibration constant for each of the segment 
averages output by SAM. The processing flow of the CALM module is described below: 

1 ) Read in the output from the segment averaging utility (run after GLA02 completes). This output 
contains segment averages (maybe 20-30 per granule) at the two calibration heights. For each 
segment average, there is maybe 10-20 latitude/longitude pairs (these are the m points along 
the orbit segment, described in 6g above). 

2) For each segment that has a valid (not invalid) PI (532), PI (1 064) or P2(532) do steps 3-6 
below. If all 3 of these are invalid, then there is no need to perform steps 3-6, below. In this 
case, we set the 3 calibration values to invalid and skip to step 9 below) 
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3) At each lat/lon point, compute the average attenuated molecular backscatter at the two 
calibration heights using ATBD equations 3.2.5 and 3.2.1 1 (here average means a vertical 
average - nominally 2 km). This requires access to the MET data at that lat/lon. 

4) At each lat/lon point, compute the ozone transmission from the top of the atmosphere to the 
calibration height (ATBD, equation 3.2.8). 

5) Compute the average attenuated molecular backscatter for the segment at the two calibration 
heights and the average ozone transmission for the segment (average of the values calculated 
in steps 3 and 4). 

6) Compute the calibration constant as the ratio of the segment signal average to the average 
attenuated molecular backscatter times the average ozone transmission (ATBD, equation 
3.2.6). 

7) Repeat steps 2-6 for each of the 20-30 segment averages. This will yield 20-30 of the 
following: Cl (532) - the lower 532 calibration constant, Cl (1064) - the 1064 calibration 
constant and C2(532) - the upper 532 calibration constant. 

8) Compute mean and standard deviation of the C values in the current granule. Throw out (set to 
invalid) those C values that are x (default x=1) standard deviations from the mean, where x is a 
variable read in from the constants file. This is done separately for each of Cl (532), Cl (1 064) 
and C2(532). Call these standard deviations cl (532). cl (1064) and c2(532). 

9) For each segment, write out to a file the following: 1 ) The start and end time for the segment, 

2) the 3 calibration values (532 upper and lower, and 1064 lower), 3) the standard deviations 
of the C values (cl (532), cl (1064) and c2(532)), 4) the three segment signal averages (532 
upper and lower, 1064 lower), 5) the segment average attenuated molecular backscatter at the 
two calibration heights, 6) the segment average ozone transmission from the top of the 
atmosphere to the calibration height, 7) the center height and thickness of the upper calibration 
zone, 8) the center height and thickness of the lower calibration zone, 9) the segment average 
532 background (B532). Note that if calibration points are thrown out during step 8 above, they 
are still output to the file, but have the value of ‘invalid’. 


The processing codes that produce GLA07 (calibrated, attenuated backscatter profiles) will then 
read the output from CALM (the calibration constants spaced roughly ‘t’ minutes apart) and will 
compute a calibration constant for each second. This process is discussed in section 3.2.11. It 
should be noted that as the laser energy decreased, the night time calculation of the calibration 
constant remained stable, but during daytime calculating the calibration constant became 
problematic for laser energy below about 5-1 0 mJ. 

The cloud search will rely on a simple threshold method, where if two consecutive bins exceed the 
threshold, then a cloud is considered found. The cloud will be output on the GLA02 product as 'the 
predicted height of first cloud top’. The cloud search is not intended to be exhaustive or the most 
sensitive. It is only meant to provide a means of detecting the first fairly dense cloud encountered. 
It will probably not be capable of detecting thin cirrus. This will be done in later processing 
(GLA09). The cloud height thus defined will be in kilometers above the local ground surface. 

The ground search begins at the end of the 1 Hz profile and works upward for a maximum of 25 
bins. The signal is searched until one bin exceeds a preset threshold value. This threshold is much 
larger than the threshold for cloud detection and was determined through simulation to be about 25 
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photons per bin. Once the ground is detected, the maximum of that bin and the following 3 is 
stored as the ‘ground return peak signal’. 

3. 1.1. 6 Saturation Flag Profiles 

The photon counting channel (532 nm) will at times become saturated by strong signals from very 
dense clouds. When this occurs, the data are no longer valid. Therefore, it is important to be able 
to recognize and flag this condition so that we can apply a correction or avoid these points in later 
data processing. The 532 channel saturation flag is in the form of a profile (SF(z)) and has a one- 
to-one correspondence with the 532 channel return signal bins. Each bin of the 532 channel will be 
checked against a maximum value (L s ) above which the signal will be considered saturated. This 
value, determined in the laboratory, is 80 counts per bin (or 156 photons per microsecond, prior to 
dead time correction). This is shown as: 

SF(z) = 0 FOR S 532 (z)<L s N sum 
SF(z) = 1 FOR S 532 (z) > L s N sum 

where N SU m is the number of shots summed for a given layer (1 , 8 or 40). Note that this test is done 
on the raw data, not the normalized signal as produced by equation 3.1.1. Parameters which will 
be read in by the algorithm and passed through as part of the output include but are not limited to: 

1 . Location of waveform peak (from altimeter channel) 

2. 1064 programmable gain amplifier setting (1 Hz) 

3. Etalon filter settings (532 channel only) 

4. Integrated 532 nm signal (P’ 532 ) from 41 to 20 km 

3 . 1.2 Error Quantification 

In this section we will try to first identify the main sources of error in the computation of normalized 
lidar signal and then attempt to quantify their magnitudes. Referring to equations 3.1 .1 and 3.1 .4, 
the main sources of error stem from incorrect knowledge of the laser energy (E) and inaccurate 
dead time correction factors for the 532 channel, and digital to analog conversion factors for the 
1064 channel. The laser energy is estimated by splitting off a small portion of the beam and 
sending it to an energy measuring device. The total energy of the beam transmitted to the 
atmosphere is then computed from this measurement. Generally, this approach to measuring the 
laser energy is accurate only to about 5 percent. The other major errors in computing the 
normalized signal for the 532 channel is the inaccuracy of the dead time correction table and the 
computation of a range dependent background for the 532 channel daytime measurements. This is 
much harder to characterize and has the added problem of changing with time. However, for 
daytime 532 channel measurements, the range dependent background is by far the largest source 
of error in the calibrated backscatter measurements. As the photon counting detectors age, and 
are exposed to continuous radiation in the space environment, their response characteristics, as 
well as the amount of detector dark current, can change. This in turn affects the dead time 
correction table. These effects are difficult to quantify, but it was observed that the 532 channel 
detector dark current did not appreciably change during the 6 years of operation. 


18 



Other factors affecting data quality are laser performance, boresite accuracy (alignment of the 
transmitted beam and the telescope field of view) and, in the case of the 532 channel during 
daytime, how well the background as a function of range has been computed and subtracted. Of 
course, the main factor affecting data quality for this mission is the laser energy loss with time. This 
applies to both the 532 and 1064 channels. We found that the boresite of the 532 channel was very 
stable and that the boresiting procedure to peak up the signal was run infrequently. The 1064 
channel did not have a boresite capability and suffered occasionally due to the laser footprint (spot 
on the ground) drifting from the telescope field of view. Likewise, signal loss will occur if the etalon 
filter is not tuned to the laser frequency. However, we found that the drift of the 532 channel etalon 
was not a problem. In the next section we will develop a set of confidence flags which are intended 
to provide a measure of data quality. 

3.1.3 Confidence Flags 

Confidence flags are meant to give an indication of data quality and our confidence that the data 
are at a level where all science objectives can be met. As mentioned above, there will be 
circumstances where the caliber of the data is reduced due to a variety of causes. Since laser 
energy proved to be the biggest factor in determining data quality, confidence flags based on laser 
energy are constructed for each channel in the LI A processing. These flags are then passed along 
and stored on the level IB and 2 products. After the laser energy drops below about 14 mJ, 
daytime measurements start to degrade. Nighttime data is still good down to about 4 mJ. 


Table 3. 532 Channel quality flags 


Laser Energy (mJ) 

Flag Value 

Comment 

>= 24 

0 

Excellent 

>= 14 <24 

1 

Good 

>= 4 < 14 

2 

Fair 

<4 

3 

Poor 


3.2 Attenuated Backscatter Cross Section (GLA07) 

3.2.1 Theoretical Description 

3.2.1.1 Overview of Processing 

The attenuated backscatter cross section falls easily from the normalized lidar signal developed in 
section 3.1 . Essentially, the only computation required to obtain the attenuated backscatter is the 
calculation of the lidar calibration constant (C). The calibration constant can be approximately 
obtained from first principles (equations 2.2 and 2.3), but, at least for the 532 channel, in practice it 
is much easier and in the long run more accurate, to obtain C from the data itself, provided 
sufficient signal is available. This approach is beneficial because it overcomes the problems 
associated with instrument drift and is self-regulating. The 532 photon counting channel has 
adequate signal for the computation of C from the data itself, but this is not possible for the 1064 
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channel. Instead we have used the results of two validation missions to calibrate the 1064 channel. 
These occurred during laser 1 (February - March, 2003) and laser 2 (October - November, 2003) 
operational periods. Utilizing the CPL onboard the NASA ER-2, spatially and temporally coincident 
measurements of clouds and aerosols were used to obtain the GLAS 1064 channel calibration 
based on the CPL’s 1064 nm measurements. The CPL 1064 channel is very accurately calibrated. 
In addition, when 532 nm was of high enough quality, we used 532 nm calibrated measurements of 
cirrus backscatter to calibrate the GLAS 1064 channel. 

Recall that what we called ‘calibration pre-processing’ (see the descriptions of SAM and CALM in 
section 3.1 .1 .3) computes the calibration constant at points ‘t’ minutes apart along the orbit. This is 
a change implemented in version 4.1 of the ATBD. The first thing which GLA07 must do (prior to 
any other processing) is to read the file output by CALM and compute the calibration constant 
which will be used for each second of the current granule. This process is implemented as a 
subroutine that is called once by GLA07 at the very beginning of its processing for a given granule. 
For ease of discussion, we will call this the calibration fitting module (CFM). The function of the 
CFM routine is simply to read in the output of CALM (C values) for the current granule and the C 
values within 1 hour prior to the start of the current granule and apply some type of fitting 
procedure to the points to obtain a C value for each second of the current granule. These values 
are then used by GLA07 to compute the attenuated, calibrated backscatter as per equations 3.2.12 
and 3.2.13. The following steps summarize the procedure: 

1 . ) Read in the output from the CALibration Module (CALM). 

2. ) Based on the value of a flag (0/1 , meaning no/yes), eliminate all C values (532 nm upper and 

lower and 1 064) with corresponding 532 nm background values greater than a threshold value 
(default = 2 photons/bin). Both the flag and the threshold value are read in from the constants 
file. 

3. ) Compute the ratio of the 532 nm calibration constant computed at the lower calibration height 

(about 1 0 km) to the 532 nm C value computed at the upper height for all points in the granule. 
If this ratio is less than 1 .0 - x, or if this ratio is greater than 1 .0 + x, then eliminate (throw out) 
the corresponding 1064 C value. Where x is read in from the constants file and has a default 
value of 0.10. Note that when the 532 nm C value at the lower calibration height is missing 
(invalid), then the value used for it in computing the ratio will be zero. 

4. ) Read in the C values that are within 1 hour prior to the start of the current granule. We now 

have x upper 532 nm Calibration points, y lower 532 calibration points and z 1 064 calibration 
points. 

Based on the value of a flag variable (0,1, 2, 3), which is read in from the constants file, the 
code then does the following: 

Flag=0: Granule mean 

5. ) Compute a mean of the x upper 532 calibration values and the z 1 064 calibration values and 

use the mean for the whole granule. A granule is two orbits. 

Flag=1: Granule mean 

6. ) Compute a mean of the y lower 532 calibration values and the z 1064 calibration values and 

use the mean for the whole granule. 

Flag=2: Running smoother 
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7. ) Using the x upper 532 calibration points, and the z 1064 calibration points compute an m point 

running average. Linearly interpolate between smoothed points to obtain the calibration 
constant for each second. The number of points to use in the running smoother (m) is read in 
from the constants file. A default value for m is 3. 

Flag=3: Running smoother 

8. ) Using the y lower 532 calibration points, and the z 1064 calibration points compute an m point 

running average. Linearly interpolate between smoothed points to obtain the calibration 
constant for each second. The number of points to use in the running smoother (m) is read in 
from the constants file. A default value for m is 3. 


In practice, to compute the 532 nm channel calibration for each second we used the running 
smoother described in (7) above. Note that there will be a default calibration value for the 1064 
channel that will be read in from the constants file. If this number is negative, then the software will 
use the calculated value of the 1064 calibration constant. If this default calibration value read in 
from the constants file is positive, then use it for the computation of 1064 calibrated backscatter, ie. 
do not use the calculated value. No such default mechanism is implemented for the 532 channel. It 
has turned out that in practice, the calculated value of the 1 064 calibration constant was never 
used. We used values obtained from co-located CPL under flights and the 532 nm measurements 
of cirrus cloud backscatter. 

Next, GLA07 computes the calibrated attenuated backscatter ((3') for both channels at 5 Hz and 40 
Hz, and correct the 532 channel (3' for times when it became saturated. Another important function 
that GLA07 will perform is the vertical alignment of the data so that each bin is referenced to height 
above mean sea level. The data acquired by GLAS (as well as the data output from GLA02) range 
in height from 41 to -1 km for the 532 channel and 20 to -1 km for the 1 064 channel. This height is 
with respect to the height above the local topography at the sub-satellite point. This is based on a 
DEM onboard the spacecraft which can have different values for each second of lidar data as 
discussed in section 3.1 .1 .3. The equations which are evaluated onboard the spacecraft each 
second to calculate the 532 nm channel (PC) and the 1064 channel (CD) range gates at which to 
start taking data are: 

PC = H sat -[H mm -H offmm ] + PC bias 
CD = H sat - [77 min — H 0 ff m [ n ] + CD bias 

where H sa t is the height of the spacecraft (from the onboard GPS which is referenced to the 
ellipsoid), Hmin is the DEM minimum, Hoffmin is the offset associated with Hmin and Pcbias and Cdbias 
are the offsets to apply to the 532 and 1064 channels, respectively. Hoffmin is set to a default of 
1 .125 km. The PC and CD biases will usually be -41 km, but can be used to move the profile either 
up (when made less than -41 km) or down (when made greater than -41 km). These will only be 
changed (from -41 km) for off-nadir pointing. The PC and CD values effectively represent the 
distance from the spacecraft to the top of the data. In Figure 3.2.1 below, this range is denoted as 
Ro. These equations are evaluated in real time aboard the spacecraft and the results are sent 
down in the telemetry data. Note that the only difference between the two equations is the bias 
term, which can be different for each channel. Also note that even though the cloud digitizer board 
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begins taking data at the same height (41 km above the local DEM value) as the photon counting 
channel (assuming PCbias = CDbias), the flight software will only send down in the telemetry those 
1064 nm data beginning 268 bins from this point (20.58 km). 


GLAS 


Let z 0 be the height above mean sea level of the top most bin of a lidar 
profile which has an arbitrary off-nadir pointing angle (0). Assume we 
know the range from the spacecraft to the top most bin (Rq), then z 0 = 
H sat - R 0 cos(0), where H sat is the spacecraft altitude in km above mean 


Let bl be the bin index of a vertical reference profile (Pj), and b2 the bin 
index of the shot at angle 0 (P 2 ). For the case where z 0 < 41 km, bl = (41 
- z 0 )/dz, where dz is the bin size (76.8 m) and b2 = 0 (ie. the top most 
bin of the profile). In the case where z 0 > 41 km, bl = 0, and b2 = (z 0 - 
41)/(dz*cos(0)). For z 0 = 41 km, bl=b2=0. 


To obtain the proper vertical alignment, we then 
proceed down the profiles bin by bin with the logic 
below to fill the vertical reference profile (Pj). 



While (bl < 548 and b2 < 548) do 

Zj — 4 1 - (bl*dz) 
z 2 = z o - (b2*dz*cos(0)) 

If ((z 2 - Zj) < dz/2) then 

Pj(bl) = P 2 (b2) ; bl++ ; b2++ 


Pj(bl) = (P 2 (b2)+P 2 (b2+1)) / 2 
bl++; b2 = b2 + 2 
endif 
endwhile 


Figure 3.2.1 The logic and algorithm that is used to correct for vertical errors introduced by off- 
nadir pointing, which may at times approach 5-1 0 degrees, and the correction to account for the 
vertical shifting of profiles from second to second as described in the text. 


Referring to Figure 3.2.1 , this means that the range from the spacecraft to the top most bin of the 
lidar profile (Ro) can potentially change from second to second, especially over areas of varying 
terrain. Thus, the same lidar bin number can correspond to different heights above mean sea level 
from second to second. The data is shifted in the vertical to account for this. GLA07 must know the 
altitude of the spacecraft and the range from the spacecraft to the top bin of the lidar profile. 
Another factor that must be considered for the GLA07 processing is the pointing angle of the laser 
beam. Normally, GLAS will be pointing very close to nadir, with pointing angles less than 1 degree. 
In this case, the effect of the pointing angle on the vertical position of the lidar return bins can be 
neglected. There will, however, be times when the pointing angle exceeds 2-3 degrees and may 
(very infrequently) be as high as 10 degrees. The effect of pointing off nadir is to cause the vertical 
distance covered by a lidar range bin to decrease by the cosine of the pointing angle. If this effect 
is neglected for larger off-nadir pointing angles, it will cause a misalignment of the bins in the 
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vertical. We will therefor take the pointing angle into account when we vertically align the data and 
place it in a coordinate frame referenced to mean sea level. 

A solution to the vertical alignment and angle correction problem is shown in the figure above. It is 
a simple algorithm and can be applied to all the data, regardless of the pointing angle. Referring to 
Figure 3.2.1 , most of the time zo will be above 41 km, as this corresponds to the situation for land 
elevations greater than zero. Over mountainous regions and high plateaus, zo can be considerably 
higher than 41 km. Over mount Everest, for example, zo will be about 50 km. In these cases, the 
portion of the lidar profile (P 2 ) above 41 km is discarded and the bottom portion of the vertical 
reference profile (Pi in the figure) is padded with a missing value. This latter point is not a concern 
as the bins to be padded are all below the ground surface (assuming of course that we have 
accurate values of A s and Ro, and that the spacecraft timing circuitry has done its job properly!) For 
the cases where zo is less than 41 km, which should not happen very often (these would 
correspond to places where the topography is below sea level, or the times of appreciable off-nadir 
pointing), the top portion of the reference profile is padded with molecular backscatter data and the 
end bins of the lidar shot (P 2 in the figure) will be discarded (as they may be below -1 km, 
depending on the pointing angle). 

3.2.1.2 Calculation of the Lidar Calibration Constant 

Since the signal return which is used in the computation of C is from purely molecular scattering, 
and the atmospheric density at these altitudes is very low, the return signal is very weak. 

Therefore, one must first integrate the return signal through a layer 2 kilometers thick centered on 
the calibration height and then average over a sufficient time span to insure adequate signal to 
noise for the computation of C. It was found that 5 to 10 minutes of averaging provided sufficient 
signal for stable 532 nm calibration calculation. The 1064 channel cannot be calibrated in this way 
and instead we use the method discussed in section 3.2.1 .1 for 1 064. Our approach for the 532 
channel is to calculate C continuously along the orbit using segments of data 1 0 minutes in length. 
The nighttime calculation is more accurate and stable (because of the lack of background signal), 
thus, it is desirable to flag each C value as being calculated during night, day or indeterminate. This 
can easily be done by looking at the background during the time that C is being calculated. Recall 
that the average background for the calibration segment was output in the calibration pre- 
processing file output by GLA02 (section 3.1 .1 .3). Based on the magnitude of the background, we 
will classify each calibration as being produced from nighttime or daytime data. If the average 
background value for a given segment is greater than about 2 photons per microsecond, then it 
can be safely assumed that it is daytime. A background less than 1 photon per microsecond would 
indicate nighttime conditions and in between would be labeled indeterminate. The values for the 
flag are -1 = night, 0 = indeterminate, and 1 = day. 


A requirement for the calculation of C is a knowledge of the average molecular backscatter cross 
section through the calibration layer and the transmission (including ozone for the 532 channel) 
from the top of the atmosphere to the calibration height. The molecular backscatter cross section 
will be needed in other GLAS processing modules in the form of profiles with the same vertical 
resolution as the lidar data (76.8 m). Thus, they will be computed in GLA07 as complete profiles 
from 41 km altitude to the surface with a 76.8 vertical resolution. This requires knowing the 
atmospheric temperature and pressure at a vertical resolution of 76.8 meters (the lidar bin size). 
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The pressure, temperature and relative humidity along the flight track are calculated from the 
ancillary MET data which is available to the GLAS ground processing system or from standard 
atmosphere tables (in the case of the 30 km calibration height). The MET data are reported at 
standard pressure levels which include temperature, relative humidity and the geopotential height. 
The geopotential height must first be converted to the equivalent geometric height and then the 
pressure (P(z)), temperature (T(z)) and relative humidity (R(z)) calculated for the bins (heights) 
between the standard pressure levels. This is accomplished with the hypsometric formula. From 
the calculated temperature and pressure profile, the molecular number density (N(z)) is calculated 
from the ideal gas law as: 

(3.2.1) N(z) = P(z) /(kT v (z)) 

where N(z) is in units of molecules per cubic centimeter, k is the Boltzmann constant for dry air in 
units of ergs per degree per molecule, P is the atmospheric pressure in units of ergs per cm 2 , and 
T v is the virtual temperature in degrees Kelvin. This equation is very similar to the equation to 
compute atmospheric density (p(z)), which is the same as 3.2.1 except that the Boltzmann 
constant is replaced by the ideal gas constant for dry air (R), which has a value of 0.0028769 m 2 S' 2 
°K' 1 . Note that the l-SIPS code must compute p(z) because it is needed for the computation of 
ozone transmission, in equation 3.2.7. The effect of moisture on atmospheric density is included 
through the use of the virtual temperature in equation 3.2.1 , but these effects are generally 
negligible above the lower troposphere. T v is computed from the relative humidity (obtained from 
the MET data) by first converting it to water vapor mixing ratio. To accomplish this, we need to first 
compute the saturation vapor pressure (e s ) which is a function of the atmospheric temperature (T) 
as: 


(3.2.2) e s = 0.61 12e 17 ' 67J7(r ~ 29 ' 66> 

and from that compute the saturation mixing ratio (q s ): 

(3.2.3) 9 S =0.622^ /(T 5 / 10.0) 


where P is the atmospheric pressure in millibars. The relative humidity is simply the actual 
atmospheric water vapor mixing ratio divided by the saturation mixing ratio times 100. Thus, the 
actual atmospheric water vapor mixing ratio is given by q = rq s / 100.0 where r is the relative 
humidity. And finally, the formula to compute the virtual temperature (T v ) is: 


(3.2.4) 


T v = 


T 

1.0-3^/5 


Following Measures (1984), from the atmospheric molecular number profile, the molecular 
backscatter cross section (p m (z,A.)) in units of mMsM is then: 

(3.2.5) p m (z, A) = 5.450/V(z)(550.0/ A) 4 10" 26 
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where X is the wavelength in nanometers (532 or 1064 nm in our case). The computation of the 
calibration constant then is: 


(3.2.6) C A =P’ JL ( Zc )/(/?Jz c ,A)T 2 (A)) 

where P' x (z c ) and P(z c ,X) 


are the horizontal average (through the calibration latitude band) of the vertically integrated 
normalized lidar signal (output from GLA02) and molecular backscatter through the 2 km thick 
calibration layer, respectively. The length of the horizontal average is defined as input to GLA02 
(default of 10 minutes). In equation 3.2.6, T 2 (X) represents the two-way path transmission from the 
top of the atmosphere to the calibration height and is composed of Rayleigh and ozone 
components as: T 2 (X) = T 2 m (X)T 2 0 (X). In this discussion, we are assuming no absorption due to 
aerosols. The ozone absorption is negligible at 1064 nm, but is large enough to consider for the 
532 channel. The ozone transmission, T 2 0 (X,z), is calculated using ozone mixing ratios obtained 
from a climatology provided by G. Labow (NASA-GSFC Code 916, unpublished data). The ozone 
mixing ratios (kg/kg) are obtained from lookup tables. The lookup tables will be grouped together 
into 10 degree latitude bands and month of year. The ozone profiles are gridded at the standard 
GLAS altitude resolution, with the first bin at 59.9796 km, stepping down by 0.0768 km to the last 
bin at number 795. 


The ozone mass mixing ratios, ro(z), are first converted to column density per kilometer (atm- 
cm/km), s o(z), using the following equation, 


3.2.7 


r o( z )p s (z) 
2.14148 x 10~ 5 


where z is the altitude in km, and p s (z) is the atmospheric density at z (obtained using the MET 
data already calculated). The next step is to calculate the ozone transmission term. T 2 0 (X,z) is 
calculated using the following equation, 


3.2.8 


T o(X,z) = exp 


-2 • c, 




s 0 (z') dz' 

glas- altitude 


where co(^) is the Chappius ozone absorption coefficient in cm- 1 . The ozone absorption coefficient 
is obtained at the correct wavelength from a table compiled in Iqbal [1984] using data from Vigroux 
[1953], co{X) is 0.065 cm- 1 at 532 nm and is zero at 1064 nm. 

The ro(z) lookup table for the 0° to 10° N latitude band is displayed in Figure 3.2.3. Ozone column 
density profiles, s o(z), were estimated for the month of July over both the equator and the south 
pole using standard density profiles [McClatcheyetal., 1971] and ro(z) from the lookup tables. The 
results are shown in Figure 3.2.4. 
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To calculate the molecular transmission, we first compute the molecular extinction profile (a m (z,X)), 
by multiplying the molecular backscatter cross section by the molecular extinction to backscatter 
ratio, which is known theoretically to be 871 / 3 . 

(3.2.9) <j m (z,A) = 8ft/3 m (z,A)/3 

The molecular optical thickness from the top of the profile (ztop) to height z is equal to the integral 
of the molecular extinction profile as shown in equation 3.2.10. 

z 

(3.2.10) T m (z,y l)= jcr m (z,A)dz 

ztop 

and finally, the two-way molecular transmission (T 2 m ) between ztop and any height z is: 

(3.2.11) r m 2 (z,y l) = e~ 2T ” (zA) 

For the atmosphere, T 2 m (z,l) is very close to one for altitudes above 1 5 km, especially at 1 064 nm 
(see Figure 3.2.5). At 9 km, the two-way molecular transmission is about 0.95 at 532 nm and 0.99 
at 1064 nm. Thus, we can assume that the two-way transmission is unity for the 1064 channel at 
the lower calibration height, but we must use the value of 0.95 for the 532 channel. Deviations from 
a purely molecular atmosphere (from aerosol above the calibration height) will lead to error in the 
assumed value of the two-way path transmission and thus to error in the calculated calibration 
constant (see section 3.2.2). 


Transmission Profiles, Molecular 



T 2 


Figure 3.2.5. The two-way molecular transmission at 532 nm (left set of curves) and 1064 nm for 
various standard atmospheres. 
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In the actual implementation of the GLAS data processing system, profiles of attenuated molecular 
backscatter (the denominator in equation 3.2.5) will be generated on a continuous basis based on 
either interpolated MET data or standard atmosphere tables which correspond to the spacecraft 
location (i.e. tropics, mid-latitude, arctic, etc). As an example, Figure 3.2.6 shows the attenuated 
molecular backscatter profiles (not including ozone absorption) for US Standard, Arctic-winter and 
Tropical atmospheres. 


Molecular Attenuated Backscatter Coefficients 



Backscatter coefficient(1/km-sr) 

Figure 3.2.6. Profiles of the attenuated molecular backscatter cross section (p m T 2 m ) at 532 nm for 
three standard atmospheres. Note that the tropical atmosphere curve is denoted by the long 
dashed curve. 


After we have computed the calibration constant at all of the points (about 8-10) along the orbit that 
were defined by the GLA02 processing, the next step is to define a calibration constant to use for 
each second. Two approaches are suggested here, but after we gain experience with the data, we 
might alter the method. For now, we will 1) calculate the average of all the calibration constants 
available for the current granule and use that one value for the entire granule or 2) linearly 
interpolate between points to obtain a unique calibration constant for each second of the granule. 
Note: the length of a granule for GLA02 and GLA07 is assumed to be two orbits. The value of a 
flag will determine which of the methods is used. Calflag = 1 means to linearly interpolate, and 
Calflag = 0 means to use the average calibration value. 
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Once the calibration constant is calculated, it must be applied to the data to obtain the calibrated, 
attenuated backscatter cross section (P’s 32 (z) and ( 3 ’io 64 (z)) for the two channels as: 


(3.2.12) B' 532 (z) = P^ 2 (z)/(C 532 T 2 o(532,z)) FOR SF(z ) = 0 

(3.2.13) j3; 064 (z) = p; 0M (z)/c l064 

Note that the ozone transmission in equation 3.2.12 is from the top of the atmosphere to height z. 
This is different than the ozone transmission used in equation 3.2.6, which is the ozone 
transmission from the top of the atmosphere to the calibration height (z c ). Note also that equation 
3.2.12 is used only if the saturation flag (SF(z)) is zero, meaning that the 532 photon counting 
channel was not saturated (as determined in GLA02). If the data were saturated, then we estimate 
the 532 backscatter from the calibrated 1064 backscatter. While this procedure can give us a 
useable estimate of the 532 backscatter, it is not entirely accurate because the magnitude of the 
scaling depends on the scattering phase function of the scattering medium which is not known. 
However, a reasonably good approximation for the 532 cross section is to simply use the 1064 
backscatter cross section as in 3.2.14. This approximation can be considered accurate to within 10 
percent for both ice and water clouds. Note that the 532 channel will be saturated most frequently 
from water clouds which tend to have larger scattering cross sections than ice (cirrus) clouds. 
Theoretical simulations indicate that the 532 channel will not saturate from most naturally occurring 
aerosol plumes, but may saturate from dense smoke from large scale (biomass burning) fires. 

(3.2.14) B' 532 (z) = B[ 064 (z) FOR SF(z)> 0 

The implicit assumption here is that we have correctly calibrated 1064 data and that multiple 
scattering (in the 1064 signal) is relatively small. The 1064 channel, with its much wider field of 
view, is much more prone to multiple scattering than the 532 channel. It is mainly the multiple 
scattering that limits the accuracy of 3.2.14. We also want to implement in the code a flag which 
dictates whether or not this substitution (3.2.14) is to take place. If the flag is set to indicate that this 
is not to occur, then the 532 channel data is not replaced, but it left intact. In practice we did not 
implement this signal replacement for saturated 532 channel bins. 

The output product for GLA07 consists of 5 Hz full profiles of P’s 32 (z) from -1 to 41 km and 40 Hz 
profiles from -1 to 10 km. The former requires averaging 8 shots from the lowest layer and the 
duplication of 8 profiles from the upper layer to form one continuous profile from 41 to -1 km. An 
example of an average 532 nm calibrated, attenuated backscatter profile is shown in Figure 3.2.7. 
For the 1 064 channel, the output will consist of 5 Hz profiles of P’io 64 (z) from -1 to 20 km and 40 
Hz profiles from -1 to 10 km, again requiring the averaging of 8 profiles from the lowest layer to 
form the entire 20 km profile. 

Output from GLA07 will include the saturation flag profiles (SF(z)) for the 532 channel output as 5 
Hz full profiles from -1 to 40 km and 40 Hz profiles from -1 to 10 km. Since the former requires 
averaging of the lowest layer, SF(z) should be set to 1 (indicating saturation) if any of the 8 shots 
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that make up the average was saturated. A detailed list of additional data output by GLA07 is listed 
in section 4.2.3. 



532 nm Backscatter Coefficient 


Figure 3.2.7 An example of the GLA07 532 nm calibrated attenuated backscatter. The green line 
is the computed molecular backscatter and the black line is actual GLAS data from October, 2003. 


3.2.2 Error Quantification 

Here we identify the major sources of error in the calculation of calibrated attenuated backscatter. 
This essentially boils down to identifying the major source and magnitude of error in the calculation 
of C. For the 532 channel, C is computed from the atmospheric scattering at specific heights 
(Equation 3.2.3). The error in 3.2.3 comes from two major sources. The first is the assumption of a 
purely molecular atmosphere in calculating the two-way transmission from the top of the 
atmosphere to the calibration height (T 2 (z c )). At the 35 km height this is ok, but the lower one goes, 
the higher the probability that some aerosol will be present. Normally, this is small since most of 
the aerosol is confined below 10 km. However, during episodic volcanic eruptions, a significant 
amount of aerosol can be injected into the lower stratosphere. Thus, the magnitude of this error will 
vary in space and time and is difficult to quantify. However, in most situations, this error will be 
negligible at the 35 km calibration height, and less than 5 percent for a 9 km calibration height. 
Further, at the lower calibration height, it will be necessary to identify and eliminate the occurrence 
of clouds in the data segment that is used to calculate C. While it is easy to find and eliminate 
dense clouds, it will be difficult to locate very thin cirrus or aerosol layers. 
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Another problem that can occur in the calculation of C is the error involved in computing the 
molecular backscatter cross section ((3m(z c ,X)) at the calibration height. For instance, if the 
temperature and pressure used to compute pm(z c ,^) were in error by 2 and 10 percent respectively 
( 4.5 °K and 1 .1 mb), then the molecular backscatter cross section would be in error by 1 0 percent. 
Thus, this error is likely to be of greater magnitude than the transmission error discussed above. A 
good way to quantify this is to plot p m for various standard atmosphere models. Figure 3.2.8 shows 
a plot of the 532 nm molecular attenuated backscatter profile for the arctic winter atmosphere (solid 
line) and the tropical atmosphere, normalized by the molecular profile for the U.S. standard 
atmosphere. This shows that at about 16 km, p m calculated from the standard atmosphere can 
differ as much as 18 percent from the p m calculated from the tropical atmosphere model. 
Essentially, this is illustrating the effect that differences in temperature and pressure have on the 
magnitude of p m . For most cases, we think that the accuracy of the MET data used to compute p m 
will limit this error to within about 5 percent. 


Ratios of Attenuated Backscatter Coefficient 

Molecular, 532 nm 



Ratio 


Figure 3.2.8. Ratio of the molecular backscatter profile computed from Arctic winter (solid) and 
tropical (dashed) atmosphere to the molecular profile computed from the U.S. standard 
atmosphere. 

As it turned out, one of the largest error sources for the attenuated, calibrated backscatter product 
was the background problem discussed in section 3.1 .1 .2. The inability to correctly calculate and 
subtract the background over highly (albedo) varying scenes causes errors in the calibration. This 
is particularly acute when transitioning through the terminator from night to day over high albedo 
areas such as Antarctica. 

3.2.3 Confidence Flags 

Confidence flags for GLA07 are based on laser energy and are defined in section 3.1 .3. 
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3.3 Particle Layer Height and Earth's Surface Height (GLA09) 

3.3.1 Theoretical Description 

3. 3. 1.1 Cloud and Aerosol Layers Height from 532 channel 

This section will present a description of the algorithms and techniques, designated GLA09, that 
will be used to find the locations of particulate layers. The GLAS 532 nm atmospheric channel 
signal will be used to locate the vertical positions of horizontal surfaces of both cloud layers and 
aerosol layers. A generalized technique will be applied initially to locate boundaries of layers of 
both particulate categories. Then, each detected layer will be determined to be cloud or aerosol 
using an objective discrimination algorithm. The cloud layer height information (up to 10 layers) is 
stored on the GLA09 product while the aerosol layers are stored on GLA08. A description of the 
algorithm to find the location of the earth's surface from the lidar return (ground signal) will also be 
given. The algorithm to detect particulate layers was originally designed to work using the 532 
channel data only. As problems arose with the 532 laser energy, it became apparent that it was 
necessary to utilize the 1064 channel for cloud detection to maximize the amount of information 
from the mission. The algorithms used to produce the cloud and aerosol layer heights from the 
1064 channel are discussed below in section 3.3.2. 

Cloud particles are those atmospheric constituents that are composed primarily of H 2 O and that are 
formed by condensation of atmospheric water vapor around condensation nuclei. Cloud particles 
form in volumes where relative humidity is at or near 100%. Cloud particles can be either liquid or 
ice and both phases can exist together. Cloud particles can condense or evaporate quickly as a 
result of the surrounding environment. Liquid droplets can exist in a super-cooled state. Clouds are 
aggregations of these particles. The aggregations typically have a layered structure, as in stratus, 
or a towering structure, as with cumulus. The two types can exist together and often a cloud has 
characteristics of both structures. A given location may be cloud free, clear, or be occupied by one 
or more types of clouds. Often, the combination of cloud types is quite complicated. Liquid water 
droplets are approximately spherical in shape. The shapes of ice particles are controlled by the 
effects of temperature, humidity, and local dynamics upon the crystalline structure. Cloud particle 
sizes usually extend over a particle size spectrum. 

We define aerosols as any particulate matter suspended in the atmosphere that is not considered 
to be cloud. The chemical composition of aerosols is quite varied. Sulfate, nitrate, salt, sand, 
smoke, and dust particles are common aerosol particles. Base aerosol particles tend to have a 
relatively long residence time compared to clouds since they are not the product of a change of 
phase of H 2 O. Base aerosol particles can be hygroscopic in which case they can serve as cloud 
particle condensation nuclei. Aerosols exist in volumes where relative humidity ranges from 0% to 
100%. In general, aerosol particle size spectrums have a smaller average size than cloud size 
spectrums. Aerosol and cloud particles coexist and often the distinction between the two is 
arbitrary. 

The distinction between cloud and aerosol layers based upon lidar backscatter signals alone is 
quite problematic. Generalized lidar signal characteristics will be assigned to each category as a 
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means to separate them in an objective scheme. It is expected that the categorization 
determination will be improved as experience is gained in analyzing space borne lidar data. 

For our purposes, we consider the layer structure to consist of a specific number of layers at any 
location. Each of these layers is a region of particles defined by a top boundary and a lower 
boundary. The lower boundary of a fog layer or a planetary boundary layer is the surface of the 
earth. A boundary exists where the density of particles exceeds an arbitrary threshold, which 
serves to define clear air. A region between top and bottom boundaries of a layer contains cloud or 
aerosol particles that could have either homogeneous or inhomogeneous characteristics. 

Because of the additive nature of scattering, cloud or aerosol atmospheric regions have greater 
volumetric backscatter coefficients than clear regions. In clear regions, radiative scattering stems 
entirely from air molecules; it is referred to as Rayleigh scattering. When particles are present, 
scattering is increased above Rayleigh scattering values. It is this enhancement in the scattering of 
photons in the lidar pulse that provides a signal that can be used to delineate particle layers in a 
lidar profile. Since absorption by water at the GLAS lidar wavelengths is negligible, the backscatter 
coefficient in particle rich regions always exceeds the Rayleigh backscatter coefficient. Because of 
this, a vertical profile of Rayleigh backscatter coefficient could be established as a baseline 
threshold to distinguish particle regions in a profile. This would be convenient since the profile can 
be readily computed when the air density is known. However, attenuation of the lidar pulse by 
intervening layers reduces the lidar backscatter signal from any given volume. Therefore, the 
Rayleigh backscatter coefficient profile can serve as only an upper limit of threshold. 

Figure 3.3.1a) provides a conceptual view of a representative lidar profile of attenuated backscatter 
coefficient together with a profile of Rayleigh backscatter. The profile was fabricated by applying 
the basic lidar equation to an arbitrarily specified atmosphere and using the GLAS lidar system 
specifications to characterize the measured signal. Cloud boundaries are clearly evident from a 
visual inspection of the lidar profile. One's perception of the profile is such that the signals above 
and below a layer provide a threshold against which the protrusion of the cloud signal is compared. 
Even where the cloud density increases gradually, such as in the cirrus layer at about 8-km, the 
boundary can be discerned to within one or two sample elements. A profile characteristic that 
masks a weak cloud boundary is the random noise superimposed upon the basic signal. The signal 
from the second layer (from the top of the profile) of cirrus is diminished because of the attenuation 
of the first. The signal from the stratus layer at 1 km is very much lessened by attenuation. Also, 
notice how the (lidar) molecular signal is diminished by attenuation in the region between 8.0 and 
10.5 km and below 6.0 km. Despite reduction of the signal due to noise and attenuation, the 
locations of cloud layers are evident. The task of an objective algorithm is to mimic what is 
perceived by eye. Figure 3.3.1 b) shows a simulated profile that has a layer with typical aerosol 
scattering characteristics between 14 km and 16 km. 

An examination of cloud signatures in lidar profiles summarized above leads us to the assertion 
that an algorithm to find cloud boundaries in lidar profiles should use localized segments of small 
signal as a baseline in testing for cloud signals. By using the profile itself, rather than a threshold 
based upon some a priori determination, we can bypass the complications that arise from the many 
different atmospheric and background conditions that will be encountered by GLAS. Also, the 
threshold can be made to be a function of altitude, which permits using values that are more 
attuned to the different types of cloud and aerosol layers at various heights. Such an algorithm can 
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be designed to be an approximation of the results that would be attained from a visual inspection of 
a profile. 

A positive attribute of an algorithm whose threshold is derived from the profile is that it can be 
implemented with very efficient computer code. The techniques required to find localized 
minimums are elementary. Only a small amount of coding is required and the solutions can be 
computed quite quickly. This will permit cloud boundaries to be found operationally at the highest 
resolution produced by the lidar. The following presents a detailed description of the algorithm. 




tfl/km-ir) 


Attenuated Backscatter Coefficient 1/km-sr) 


Figure 3.3.1. a) Simulated GLAS profile in a cloudy atmosphere. Two cirrus layers and one stratus are present. 
The optical depths are from top to bottom 0.5, 1.0, and 1.5; b) simulated profile that includes a layer with 
stratospheric aerosol scattering optical properties between 14 km and 16 km. 


Cloud and aerosol layer boundaries will be found at four time resolutions. These are, from coarsest 
to finest, 0.25 Hz, 1 Hz, 5 Hz, and 40 Hz. To do this, the GLAS time series will be divided into a 
sequence of independent 4-second segments. These segments will be subdivided into four 1 
second segments. Each of these will be divided into 5 segments and these will be divided into 8 
segments, which will occur at the basic GLAS 40 Hz. frequency. Profiles of attenuated backscatter 
coefficients will be produced at 40 Hz and 5 Hz by GLA07 and serve as input into the cloud 
boundary algorithm. The 1 Hz and 0.25 Hz profiles will be produced by averaging the higher 
frequency data. 

Boundary search operations will be applied to 0.25 Hz profiles first. Results at finer resolutions will 
be made only in vertical regions where layers were detected at a coarser resolution first. The 
reason for this procedure is that the smaller signal to noise characteristic at higher resolutions will 
tend to obscure any layers not detected at lower resolutions. This technique will fail to detect some 
cloud layers that are composed of horizontally sparse and rarefied patches. But such layers are 
presumed to be insignificant for climatological studies. 

The basic layer boundary search technique will be the same for each of the four resolutions. Since 
the 0.25 Hz resolution profiles will be those first searched for the presence of cloud layers, we will 
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focus first on those in our description of the search algorithm. The finer resolutions will use the 
results of coarse resolution searches to eliminate portions found to be cloud free. 

Four one second attenuated backscatter coefficient profiles will be averaged together to produce a 
four-second averaged profile. A discussion of the potential difficulty caused by varying ground 
height among the four one second profiles will be given in a later section. The profile will be divided 
into a small number of segments. The optimum number will be found by applying the technique to 
simulated and proxy data sets to determine the means to obtain the best results. The number will 
likely be in the range of five to ten. The objective is that each segment has some samples that are 
in particle free portions of the profiles. A characteristic signal from particle free segments can 
reasonably serve as a layer signal threshold. In general, it will not be known, a priori, whether a 
segment has layer free samples. The difficulty is that rarefied layers are not easily discerned in a 
noisy profile. Each of the segments will be searched for its minimum value. Also, in order to better 
characterize each segment, the mean and variance of the sample values will be computed for 
each. In the cases where a segment has particle free regions, the minimum values will represent 
the attenuated signal from atmospheric molecules with negative random noise excursions 
superimposed. These will thus represent the absolute minimum that any layer-distinguishing 
threshold could be in each of the segments. A reasonable maximum threshold would be the 
computed molecular backscatter coefficient. Together, these values represent a range of values 
that could serve as layer signal threshold. To employ the molecular backscatter coefficient profile 
as the maximum threshhold value, a truncated lidar signal profile will generated. This profile will 
have the molecular profile as the upper limit of its value at any height. 

To find an optimum threshold value within the threshold envelope, it is necessary to find a measure 
of random noise because the lower limit boundary of threshold values is strongly influenced by the 
magnitude of random noise. This magnitude can be represented by the standard deviation of the 
lidar signal in a layer-free profile segment. Based upon our experience, we can assert that the 
atmosphere is, in general, free from non-molecular, strong-scattering species in the 18-19 km 
layer. Therefore, the noise of the lidar signal there stems mostly from the molecular scattering 
signal and the background energy. The variance of the truncated lidar signal from this high region 
will be used as a measure of the random noise contained in the lidar signal. 

Once a typical molecular signal variance has been computed, layer signal thresholds can be 
computed for each of the profile segments. In each segment, the threshold will be the sum of the 
minimum and a constant fraction of the square root of the variance. In some cases, the sum would 
exceed the computed molecular signal The value of the fraction will be determined from GLAS 
signal modeling studies but it will likely have a value in the range of 0.25-0.5. A profile of layer 
signal threshold will be then constructed by piecemeal, linear interpolation of the segment values. 
The interpolation would be done at GLAS vertical resolution. The interpolated profile will serve as a 
layer signal baseline upon which the presence of layer signals will be tested. 

The threshold profile described above will have the following positive attributes: 1) threshold values 
will be computed from the profile itself and will automatically adjust to the current situation; 2) the 
threshold computed at given level will be influenced by the attenuation of the lidar signal by higher 
layers; 3) the technique will be valid for any time resolution. A negative attribute is that the 
statistical nature of the computation of variance introduces some uncertainty into any particular 
result. 
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Once the profile of layer signal thresholds is established for a lidar signal, the layer boundaries are 
sought in the following manner. Starting at the top of the profile, the lidar profile is tested on a 
sample by sample basis. If a value is found to exceed the threshold, it is deemed a potential layer 
sample. If a specified number of potential consecutive layer samples are found, the segment is 
designated a layered region. The top of the layer is located at the height where the highest of the 
consecutive samples was found. The high-to-low testing continues under the stipulation that the 
profile is in a layered segment. The layer designation continues until several consecutive samples 
are found to be less than the layer threshold. In that situation, the profile is considered to be in a 
layer free region. The bottom of the layer is the point where the first of the consecutive particle-free 
values was found. The testing continues downward for the top of another layer. The profile will be 
so analyzed for layers to the DEM based location of the earth's surface. 

The layer boundary analysis for a 0.25 Hz profile will be used as the basis for the equivalent 
analysis of the four 1 Hz profiles that it encompasses. The layers at which 1 Hz layer boundaries 
will be produced will be limited to those vertical intervals where layers are detected at 0.25 Hz. The 
reason for this design is that averaging to produce 0.25Hz profiles will result in samples with a 
large signal to noise characteristic, which will make it least likely to result in the fewest cases of 
incorrectly identifying layers. The 1 Hz data will have a smaller signal to noise ratio value. Limiting 
the results of the 1 Hz search to the layers as 0.25 Hz will minimize false layer results at 1 Hz. For 
practical reasons, the search for layers at 1 Hz will use entire 0-20 km profiles, but the layered 
regions found will be limited to those found at 0.25 Hz. The implication of these limitations is that 
any layers which are not substantial enough to produce a detectable signal at 0.25Hz are not 
considered to be significant at finer resolutions. 

The results of the search for layers from 5 Hz. profiles will be limited by the results from the 1 Hz 
profiles in a manner equivalent to the limitations imposed upon 1 Hz by 0.25 Hz. The same search 
algorithm will be applied from 0-20 km but the resulting detected layers will have to be among the 
layers detected at 1 Hz. or they will be discarded in the output. The situation for 40 H. will be 
slightly different to accommodate to relatively small signal to noise at that frequency. Layer 
detection at 40 Hz will be limited to regions where one or more layers were detected in the 5 Hz 
profiles. If one or more layers are found in a 40 Hz profile, only the lowest one will be recorded. 
This procedure will allow detection of low cloud layers that typically have strong lidar signals and 
that have horizontal distributions that vary at relatively high frequencies. 

There are difficulties that arise from the variable ground height that may exist along the distance 
interval over which the average profiles will be produced. GLAS will produce vertical profiles that 
will use the local DEM value as the reference and lower boundary. The DEM values will be 
updated every 1 second and so four DEM values will be used in the construction of the 20, 5 Hz 
profiles which will be used to produce a 0.25 Hz profile. For purposes of layer boundary detection, 
the value of the highest DEM boundary used within the 4-second interval will be considered the 
lowest altitude at which to search the profile for layers. Also, since the one-second period of the 
DEM updates will probably not be synchronized with the 1 Hz lidar profiles, the higher of the two 
DEM values spanned by the duration of the profile will be used as the lower boundary for the 
search. Individual 5Hz and 40 Hz profiles will be contained within a single DEM interval, so this 
overlap problem will not exist. 
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A very important characteristic of downward looking lidar must be noted. As the laser pulse travels 
through the atmosphere, the scattering processes diminish its energy. In the case of a relatively 
small cumulative optical depth, the reflection of the pulse from the earth's surface has enough 
energy to be detected. If the cumulative optical thickness of the scatterers is large enough, the lidar 
signal will be reduced to the background level regardless of the magnitude of backscatter 
coefficients and no ground signal will be detected. No bottom boundary can be detected. Thus, 
when no return signal is detected from the earth's surface, the height of the bottom of the lowest 
layer is, in general, an invalid value with no relationship to the actual location. If a ground signal is 
detected, the uncertainty in the location of the bottom of the lowest layer increases as the ground 
signal's strength decreases. 

3.3.1.2 Objective Layer Discrimination Procedure 

The cloud and aerosol layer boundary detection algorithm will find top and bottoms of layers of 
both categories of particles. In order to assign a type to a layer, it will be necessary to test it to 
determine how well it matches characteristics ascribed to each category. A scoring of the testing 
will indicate the likelihood that a layer is either cloud or aerosol. It will be possible to categorize 
most layers with much confidence. However, some layers will have characteristics of both 
categories and the classification of these will be much less certain. In this section, we summarize 
the characteristics of each category and indicate the tests that will be done to discriminate between 
cloud and aerosol layers. 

Table 4 lists layer attributes that can be detected with a lidar signal. Under the particle-type 
categories of cloud and aerosol is an indication of the relative values expected for each category. 


Table 4. Comparative attributes of aerosol and cloud layers. 



Aerosol 

Cloud 

Signal Magnitude(top) 

Smaller 

Larger 

Signal Gradient(top) 

Smaller 

Larger 

Altitude(top) 

Lower 

Higher 

Horizontal Extent 

Wide spread 

More localized 

Horizontal Uniformity 

More uniform 

Less uniform 

Vertical Extent 

Larger 

Smaller 

Vertical Uniformity 

More uniform 

Less uniform 

Relative Humidity 

Lower 

Higher 

Attenuation 

Lower 

Higher 


Each attribute is listed and briefly discussed below. 

a) Signal magnitude: The average and maximum signals in an aerosol layer tend to be less in an 
aerosol layers compared to cloud layers Cloud layers that have signal magnitudes of the same 
order as aerosol layers are typically cirrus layers that will be much higher in the troposphere 
than aerosol layers. 

b) Signal gradient at layer top: The top boundary of aerosol layers will tend to be less distinct, 
than tops of cloud layers. 

c) Altitude: Cloud layers with small lidar signal characteristics will tend to be higher in the 
troposphere than aerosol layers 


37 
































d) Horizontal extent: Elevated aerosol layers will tend cover larger areas. This is a weak 
distinguishing characteristic. 

e) Horizontal uniformity: Elevated aerosol layers will tend to be more well mixed and therefore 
more uniform that cloud layers. 

f) Vertical extent: Elevated aerosol layers will tend to have a larger detectable vertical extent 
than cloud layers at the altitude at which they are found. 

g) Vertical uniformity: Elevated aerosol layers will be more uniform vertically. 

h) Relative humidity: Measured relative humidities approaching 100% are necessary for the 
presence of a cloud layer. Elevated aerosol layers can exist at much lower relative humidities. 

i) Attenuation: The vertical region where elevated aerosol layers exist is the lower troposphere. 
The clouds in this region generally have a volumetric backscatter cross section that is much 
larger than aerosol layers. Consequently, the optical depth of cloud layers will tend to be much 
higher than aerosol layers. 

Table 5. Computed discrimination criteria parameters for the aerosol and cirrus layers shown in 

Figures 3.3.2 a) and b). 



Aerosol 

Cloud 

Signal Magnitude(top) 

2.1 x 10' b /m-sr 

r 3.0 X 10' b /m-sr 

Signal Gradient(top) 

-7.5x1 0'/m-sr/km 

-1.5x10' b /m-sr/km 

Altitude(top) 

-5 km 

-11 km 

Horizontal Extent 

N/A 

N/A 

Horizontal Homogenity 

0.25-0.35 

0.2-1 .0 

Vertical Extent 

-4 km 

-6 km 

Vertical Homogenity 

0.12-0.30 

0.4.0. 8 

Relative Humidity 

-35% 

>75% 

Attenuation 

0.3 

0.6 


The actual values used for the discriminating criteria presented above will be determined from 
modeling studies and from studies of atmospheric lidar data taken by the NASA ER-2/CLS and 
other high altitude and ground based lidars. In general, low level layers that display weak and 
uniform lidar signal strength characteristics in a low relative humidity environment will be classified 
as aerosol layers. Most other layers will be considered cloud layers. 

An example of the contrasting characteristics of cloud and aerosol layers is depicted in Figs. 3.11a 
and 3.11b. Table 5 shows results of rough computations of the discriminating criteria discussed 
above. 
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Figure 3.3.2 a) Satellite lidar signal Figure 3.3.2 b) Satellite lidar signal 

cross-section showing a representative cross-section showing a representative 

elevated aerosol layer. cirrus layer. 

In an operational environment, difficulties in quantifying horizontal extent, horizontal homogeneity, 
vertical extent, vertical homogeneity, relative humidity, and attenuation probably precludes using 
these characteristics to distinguish between cloud and aerosol layers. Therefore, only the signal 
magnitude, signal gradient, and altitude of the top of each layer will be used in the layer 
discrimination procedure. The following description gives the details of the discrimination 
technique. 

The discrimination algorithm is based upon a thresholding process where the value of a single 
parameter serves to distinguish between the two categories of scatterers. For any parameter, there 
will likely exist a range of values that could indicate either cloud or aerosol. Therefore, it is 
reasonable to use a measured value of the parameter to find the probability that the layer belongs 
one of the categories, for instance, clouds. The arbitrary value of the parameter also determines 
the probability that a layer belonging to the other category, aerosols, is falsely assigned to the 
cloud category. The optimum value for the parameter is that which maximizes the probability that a 
cloud layer will be correctly identified while minimizing the probability that an aerosol layer is 
identified as a cloud. In this procedure, a correct selection of a layer as a cloud is considered a true 
positive and an incorrect selection of an aerosol layer as a cloud is a false positive. 

The discrimination algorithm will be implemented in the following manner. The cloud detection 
algorithm will be applied to 0.25Hz data. Each detected layer will be assigned to an altitude 
category based upon the height of the top of the layer. For each layer, a parameter composed of 
the product of the layer's maximum signal and maximum vertical gradient magnitude will be 
computed. This product will serve as a discriminator for cloud and aerosol layers. Cloud layers will 
tend to have a significantly higher value. The value of the product will be compared to a threshold 
value previously determined for each altitude category. If the product exceeds the threshold, then 
the layer will be deemed a cloud layer. Otherwise, the layer is deemed an aerosol layer. The value 
of the threshold is arbitrary. It will be set at the lowest point where the probability of true positive of 
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cloud designations is considered high enough when balanced against the probability of false 
designation of an aerosol as a cloud. All of the layers at higher frequencies associated with a layer 
designated cloud or aerosol will be considered to belong to the same category. Therefore, the 
cloud-aerosol discrimination will need to be applied only to 0.25 Hz data. The altitude category and 
discriminator threshold values will be stored in a table that will be read when the program is 
initialized. The values will be determined from statistical studies of existing ground-based and 
airborne-based lidar system databases. All layers that are classified as aerosol are stored on the 
GLA08 product. The following diagram depicts the logical flow of the algorithm. 



Layer discrimination algorithm 
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3 . 3 . 1.3 Corrections for false positive and certain false negative results 

Random noise fluctuations will result in layers being detected where no clouds and aerosol layers 
actually exist. False positive results of this type generally have the three distinguishing 
characteristics: 1) they have a very small vertical extent; 2) they have small horizontal extent; 3) 
they have a small total integrated signal. In order to minimize false positive results, the following 
tests will be incorporated into the algorithms at each testing resolution. The vertical depth of a layer 
will be compared to a minimum value. If the thickness is less than the minimum, the integrated 
vertical signal will be computed and compared to a minimum. If the integrated signal falls below the 
minimum, the result will be considered negative. For 1Hz, 5 Hz, and 40 Hz resolutions, a test will 
be made for the existence of the layer in adjacent profiles. If the thin layer exists in less than a 
minimum number of profiles, the result will be considered negative. 

It is possible for the threshold method to fail to detect vertically thin but optically dense layers. 
These types of layers would result in a large lidar signal in a number of vertical samples less than 
the minimum necessary to be classified a layer. In such cases, no ground signal would be 
detected. Stratocumulus clouds could give such a signal. To detect such layers, an additional 
search of the signal profile will be made each time no ground signal is found. Starting at just above 
DEM level, the signal profile will be tested for large values. When a value is found to exceed a 
certain minimum, which will be much larger than the threshhold values, a cloud layer, which has a 
vertical thickness less than the threshhold algorithm minimum, will be put at that height 

The values of the parameters used in these tests will be determined from modeling studies and 
from actual data as it becomes available. 

3. 3. 1.4 Remedy for Day/Night Bias 

Reflected solar energy is the source of two major components of total lidar signals from sunlit 
regions. These are constant offset signals, which are usually referred to as background, and 
random noise fluctuations, which are measured by the square root of the variance (root mean 
square, RMS) of random noise superposed upon the profile. Both components increase as the 
strength of reflected energy increases. The background component of a GLAS signal profile will be 
determined by averaging the signal in the portion of the profile where no laser signal is present (the 
background region of a profile). The background signal will be subtracted from the total to leave 
only the laser signal and random noise to comprise the total signal. 

Our methodology to determine layer boundaries is based upon constructing a layer signal 
threshold profile where the value of the threshold is strongly dependent on the RMS value of signal 
random noise. A larger RMS value will lead to larger threshold values. As indicated above, the 
magnitude of the RMS noise will be larger, in general, during daylight observations than those 
taken in darkness. The resultant threshold values become larger. This results in the layer detection 
technique being less sensitive to a given, small layer signal during daylight observations than 
during night observations. Layers with a certain level of weak signal will be detected in night 
observations but not in day observations. A day-night layer detection bias is the result of this 
procedure. Such a bias would hamper certain types of layer studies. 
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A solution to the day-night bias is to determine a threshold profile that is diurnally invariant and use 
this profile for all layer detection operations. A constant threshold profile would eliminate the 
differences caused by changing RMS magnitude of random noise. But, in order to eliminate false 
layer detection during daylight observations, such a threshold profile would have values that are 
greater than necessary for dark observations. For nighttime application, the method would be less 
sensitive than what is possible. Significant cloud and aerosol layers that could be resolved would 
go undetected. 

In order to give both complete and unbiased layer boundary results, the GLAS algorithm will be 
applied twice. One application will use a threshold profile based upon the observed RMS noise of 
the backscatter profile (as discussed in section 3.3. 1.1). The second application of the algorithm 
will use a threshold profile based upon a diurnally invariant threshold profile. The procedure is as 
follows. The boundary algorithm will be applied exactly as described in prior sections. This 
algorithm employs a threshold profile that uses the RMS magnitude of the profile noise as one of 
its components. Detection of cloud and aerosol layers in this manner will be the most sensitive for 
a given situation. Layer locations will be found and recorded at each of the temporal resolutions 
(0.25Hz to 40 Hz). After this operation is completed, the algorithm will be reapplied, this time using 
a threshold profile that incorporates an invariant noise component. The lidar signal will be 
compared to the threshold only in portions of the profile where layers were detected using the 
variable threshold profile. If the presence of a layer is indicated during this testing, it will be 
recorded in a true/false variable but its top and bottom boundaries will not be re-computed. This 
application will proceed through each of the resolutions. The result of the dual application of the 
layer boundary algorithm will be: a) a set of layer boundaries at each of the temporal resolutions, 
determined with the variable threshold profile; b) a set of corresponding true/false flags indicating 
whether each of the layers was detected using the diurnally invariant threshold profile. 

Determination of the invariant RMS noise component will require appropriate GLAS simulation 
studies. A threshold profile must yield results where few significant layers are missed and where 
few false positive results occur. A trade-off between these two competing requirements always 
exists in finding a threshold. Modeling studies will permit the final determination of the threshold to 
be based upon the expected performance of the GLAS lidar and will permit an estimate to be made 
of the sensitivity and tolerance of the algorithm. 

3. 3. 1.5 Polar Stratospheric Clouds (PSCs) 

Polar stratospheric clouds are layers of particles that occur in Polar regions during winter seasons 
at the respective poles. These layers reside in the stratosphere from 15 to 30 km in altitude. The 
layers are composed of particles of various chemical compositions. These layers are more properly 
classified as aerosol layers than as H 2 O cloud layers. They can reside above the cloud and aerosol 
layer boundary algorithm upper limit (20km). Any PSC found as part of the layer detection 
algorithm will be classified an aerosol layer. They will be analyzed as part of the aerosol detection 
algorithm (see section 3.4. 1.2). 

3. 3. 1.6 Bottom of Lowest Layer 

A short discussion concerning the ambiguity in the altitude of the bottom of the lowest detected 
cloud layer is given in the final paragraph in section 3.3.1 .1 . Two additional assertions can be 
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made concerning this. First, if the ground signal is not detected, the bottom of the lowest detected 
layer is not determinable and additional layers may exist below the last layer. Second, the 
uncertainty in the location of the tops and bottoms of each detected layer increases as the 
cumulative optical thickness from the spacecraft increases. There is a flag on the GLA09 data 
product for each cloud layer detected which indicates whether the bottom found is considered to be 
the true cloud bottom or not. If the latter, then it means that inspection of the signal below this point 
revealed no further layers or a ground return were found. 

3. 3. 1.7 Earth's Surface Height 

The detection of the earth's surface (GLAS ground signal) presents a problem very similar to that 
of detection of layer boundaries. In fact, the algorithm is simplified because only one surface is to 
be found. Also, because the timing of the GLAS laser is synchronized with a 1 degree DEM of the 
earth's surface, the algorithm will have an approximate location available and the search can be 
limited to a small interval surrounding that height. 

The characteristics of the ground signal in a GLAS profile are affected by the time resolution of the 
profile. Since the profile samples are much larger than the length of the laser pulse, the ground 
signal will be contained in only one or two samples at 40Hz resolution. However, the effective 
ground signal can broaden when more than one laser pulse is used to generate a profile. This 
broadening is caused by the variability of ground location over the horizontal extent that is used to 
generate the profile. If the terrain is rugged, the broadening would extend over 10 or more pixels 
for a 0.25 Hz profile, which would lead to a significant ambiguity in the meaning of ground location. 
Thus, a modified definition of ground signal is required of low resolution profiles. 

Random noise can mask the ground signal. This is especially true for higher frequency profiles 
where signal attenuation reduces the pulse strength. This effect is generally less important when 
multiple shots are used to produce a profile. 

The competition between higher precision results from high frequency profiles and higher reliability 
from lower frequency profiles leads to compromise algorithm design where the 5 Hz profiles will be 
used as the primary ground-location analysis frequency. The 5 Hz results will be averaged to 
produce ground locations at 1 Hz and 0.25 Hz. In addition, the location of the 40 Hz ground signal 
will be limited to an elevation interval close to that found for the encompassing 5 Hz profile. 

The search for ground signal in a 5 Hz profile will proceed as follows. Since the GLAS laser is 
timed so that the final 1 3 samples of a profile occur after the level of the DEM elevation, the initial 
guess for the height of the earth's surface is at the 1 3 th sample from the end of the profile. In such a 
case, the signal in the final 12-13 samples would be purely background with random noise 
superimposed. This permits a ground signal threshold to be computed from the signal in this 
segment. To do this, the mean, median, maximum, minimum, and variance of the final 20 samples 
will be computed. A threshold will be computed by adding the median and the square root of the 
variance multiplied by a factor that is a function of the current conditions. The value of the factor 
will be determined from simulation and proxy-data studies, which will reveal the optimum value to 
use in different circumstances. The values of the samples, beginning with the latest and 
proceeding to the earliest (bottom to top), will then be compared to the threshold. If a single value 
or several non-consecutive values exceed the threshold by a relatively large amount (perhaps 
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three standard deviations for instance) then the earliest (lowest height) of these will be considered 
the ground signal. Otherwise, if there are one or more occurrences of one or two-only consecutive 
samples that exceed the threshold, then the lowest of these will be considered the ground signal. 
The higher sample of any ground signal pair will be selected as the ground signal. If no such 
results are found, then the ground signal will be considered undetectable for the profile. Once all of 
the 5 Hz ground signals within a 1 Hz or 0.25 Hz averaging segment are found, the detected 
ground signal heights of the 5Hz results will be averaged to produce the ground height for each of 
the lower frequencies. 

Finally, this same ground signal detection algorithm will be applied to each of the 40 Hz profiles. 
The parameters that are derived from modeling studies will have different values than those for 5 
Hz. The low signal to noise will result in a higher rate of falsely detecting ground signal. 

3.3.2 Cloud and Aerosol Layer Height from 1064 Channel 

3.3.2.1 Overview 

The 1064 cloud and aerosol detection algorithm uses the calibrated, attenuated backscatter 
profiles as computed in section 3.2 and stored on the GLA07 product. The cloud heights are 
determined at 0.25, 1 and 40 Hz. Aerosol layers are given only at the 0.25 Hz resolution and the 40 
Hz cloud layers report a cloud top, not a cloud bottom. As with the 40 Hz 532 channel derived 
cloud heights, the 1064 channel 40 Hz heights only cover the altitude range 0 to 10 km, as 40 Hz 
data is down-linked only for that height range. The 1 and 0.25 Hz cloud layer heights extend from 
the ground up to 20 km altitude since the 1064 data above 20 km is not downlinked from the 
satellite. 

The 1064 layer detection algorithm first averages 20 seconds of data and searches for the layers 
from 20 km altitude to 250 m above the surface. If a layer is found at the 20 second resolution, five 
four second averages from the 20 second interval are formed and the 4 and 1 second layer 
detection algorithm (detailed below) is applied to obtain layer heights at the 4 second resolution 
except that the search is narrowed to a band defined by the top and bottom of the layers detected 
at 20 second resolution. If layers are detected within any of the five 4 second average profiles, then 
four one second average profiles are constructed and searched for layers again limiting the search 
to a band surrounding the layers detected at the 4 second resolution (using the same algorithm). 
The width of the layer search band is defined as 500 m above and below the determination of the 
cloud top and bottom, respectively, from the coarser resolution before it. 

3.3.2.2 1064 Layer Detection Algorithm Structure 

3.3.2.2.1 20 Second Layer Detection 

For the 4 and 1 second resolutions the same 1 064 layer detection algorithm is employed. For the 
20 second average search, a separate algorithm is used and described here. The first step is to 
apply a 3 point binomial smoother to the entire 20 km backscatter profile. A layer top theshhold (T t ) 
is then computed from the sum of the magnitude of the 1064 molecular (Rayleigh) scattering at the 
current height ((3 m ) and a baseline (Y) multiplied by a factor (F) that depends on the current height 
as: 
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(3.3.1) T t = (p m (z) + Y)*F 

Where Y = 1.0x1 O' 6 and F = (0.70d0 +(i_bin_cnt/2)/180.0d0). Here, i_bin_cnt is the return bin 
number starting at 0 at 20 km altitude and reaching 260 at the ground. The signal is then tested 
starting at 19 km (i_bin_cnt = 14) altitude moving downwards for the occurrence of 3 consecutive 
bins for which the signal exceeds the top threshold (Tt). Once this condition occurs, the layer top is 
defined as the height of the first of the three bins to exceed the threshold. After encountering a 
layer and while within a layer, the bottom threshold (Tb) is computed bin by bin as we move 
downward through the layer as: 

(3.3.2) Tb = Tt — (i_bin_cnt - i_bin_top) * d_top_thr * 0.0010) 

Where i_bin_top is the bin representing the top of the current layer. Thus, Tb is initially equal to Tt 
at the top of the layer and decreases by about 1 % for every 750 m as we travel downward within 
the layer. This continues until we find 3 consecutive bins below Tb at which time the first of these is 
deemed the layer bottom. The search then continues with a new T defined from equation 3.3.1 until 
either 10 layers have been found or we reach 250 meters above the ground. If we are in a layer 
and reach the 250 m level and have not yet determined a layer bottom, the bottom height is set to 
250 meters above the local surface height determined from a 1x1 degree DEM. 

3.32.2.2 4 and 1 Second Layer Detection 

Once a layer is found at the 20 second resolution, then five, 4 second average profiles are formed 
and searched for the presence of the same layer. Here, the whole profile is not searched, but 
rather the search is limited to a band defined by the top plus 1500 m and bottom minus 500 m of 
the layer obtained from the 20 second layer search. Beginning at the top of this band and working 
downward to no lower than the height of the maximum signal within the band, a top threshold is 
formed bin by bin as: 

T, = (p m (z) + Y) * (0.90 + i_bin_cnt/200.0) 

Where Y = 1 .2x1 O' 6 and i_bin_cnt is the bin number at the current height of the search (varies 
between 1 (20 km) and 260 (about 250 m above the ground)). When 3 consecutive bins are 
greater than T t then the top of the layer is the height of the first of the 3. If the top is found, a search 
is initiated for the layer bottom beginning at the height of the maximum signal within the band 
(which is also the bottom limit for the search for the top) continuing to 500 m below the bottom 
obtained from the 20 second search. The bottom is defined as the height of the first of 3 
consecutive bins that are less than Tb defined as: 

T b = (p m (z) + Y) * (0.90 + i_b i n_to p/200.0) 

The profiles that contain layers detected at the 4 second resolution are then searched at a 1 
second resolution using the same algorithm, this time using the 4 second layer heights to define 
the band within which to search. An example of the 4 second 1064 cloud heights is shown in 
Figure 3.3.6. 

3.3.2.2.2.1 Confidence Flag for 1 and 4 Second 1064 Cloud Top 
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The confidence or quality flag for the 1 and 4 second cloud layers are as follows: value 15 = cloud 
layers were not searched for; value 0 = cloud layers were searched for but not detected; values 
1-14 indicate increasing confidence of good cloud retrieval (value 1 = least confidence, value 14 = 
greatest confidence). The values of 1 through 14 are computed from the average signal in the 
first 300 m of the cloud divided by the average signal above the cloud top. 

3.3.2.2.22 Cloud/Aerosol Discrimination for 1 and 4 Second 1064 Layers 

A rudimentary cloud/aerosol discrimination is applied to the layers found by the 4 and 1 second 
1064 searches. If a layer is deemed cloud, it is stored on the GLA09 product and if it is aerosol, it is 
stored on the GLA08 product. If a layer has been detected and the current latitude is between 60 N 
and 60 S, and the layer top is below 6 km altitude, then the following test is applied. Outside of this 
latitude altitude range, the layer is assumed to be a cloud. 

Three thresholds are defined based on the top height of the layer (z) and the layer number (N): 


Ti = 6.00x1 0-e nr 1 sM / N for 6.0 >= z >= 4.0 

T 2 = 7.50x1 0- 6 m- 1 sr 1 / N for 4.0 > z >= 2.0 

T 3 = 9.00x1 0-6 nr 1 sr 1 / N for 2.0 > z > 0.0 

If the average backscatter within the layer is less than the thresholds defined above (for the 

respective height of the layer top), then the layer is assumed to be aerosol. 

3.3.22.3 40 Hz 1064 Cloud Detection 

The 40 Hz 1064 nm cloud detection was put in place after it became apparent that the altimetry 
community would like to have information on the presence of clouds on a shot by shot basis. 

Ideally, this measurement would be made with the 532 nm channel, but because of the laser 
issues with that channel, we decided to use 1064 nm so that many more of the operational periods 
would have viable results. Even so, as the decline of the 1064 laser energy became substantial, 
this algorithm will no longer provide stable results. After the laser energy falls below about 35 mJ, 
the 40 Hz cloud detection begins to be unreliable. However, prior to that time, the 40 Hz 1064 
cloud detection is quite good. An example of the output of the algorithm is shown in Figure 3.3.7. 

The algorithm first applies a vertical and horizontal smoother to the 40 Hz 1064 profiles using a 
triangular smoother of width 3. The resulting profiles are then vertically smoothed using a box car 
smoother of width 7 bins. 3 point running vertical smoother to the 1064 nm profile and then applies 
a 5 profile running smoother. The resulting profile is searched from the top (10 km) downward to 
150 m above the DEM indicated ground. If 4 consecutive bins exceed the threshold (1.1x1 O' 5 ), then 
the first of the 4 bins is the 40 Hz cloud top. 

For all cloud layers detected, a quality or confidence flag is then computed from the average of the 
top 4 bins of the cloud layer divided by the cloud top threshold (l.lxlO- 5 ). If a cloud is not detected, 
the integrated backscatter is calculated from the smoothed profile. If the integrated backscatter is 
greater than (0.5x1 O' 5 ), then the quality flag is set to 13 and the cloud top height is set to 100 m. 

3.322.3.1 Confidence Flags for 40 Hz 1064 Cloud Top 
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Value 15 = No clouds. 

Value 14 = Indicates the likely presence of low clouds (< 150 m) based on elevated signal from the 
two bins above the ground return bin that were not detected directly by the cloud search algorithm. 
When this occurs, the 40 Hz cloud top height (i_Frir_cldtop) is set to a value 
of 0.10 km. 

Value 13 = Indicates the possible presence of a cloud based on the value of the integrated signal 
parameter (LFRirJntsig) that was not detected directly by the cloud search algorithm. When this 
occurs, the 40 Hz cloud top height (i_Frir_cldtop) is set to a value of 10.0 km. 

Value 0 - 12 = Cloud detected by cloud search algorithm with higher numbers indicating a stronger 
average signal from the region starting at cloud top and extending 500 m below cloud top height. 

3.3.3 Error Quantification 

Multiple scattering is a potential source of large error in determining the boundaries of layers and 
the earth's surface from a space-borne lidar. The multiple-scattering process causes secondary 
photons to take deviated paths back to the lidar receiver where they are combined with the single- 
scattered signal of later samples. This causes the later sample to appear to have a larger signal 
than that based solely upon the density of the scatterers. A possible result of this is that a layer's 
lower boundary is analyzed to be at a lower altitude than it actually is. Fortunately, the vertical 
resolution of the boundary analysis is, at best, 76 .8 m. Our experience with spaceborne lidar 
indicates that the multiple scattering effect is significant, at this resolution, only in dense low clouds. 
Since these clouds usually fully extinguish the laser pulse, no ground signal would be detected and 
the lower boundaries of these clouds would be unknown. Because of this, it is not expected that 
multiple scattering will have a significant effect on the quality of the results of the boundary 
algorithm for most clouds. 

The quality of the results of a layer boundary algorithm can be divided into two components: a) true 
or false determination of the existence of layers; b) precisely locating the top and bottom of layers. 
Errors in component a), designated false positive or false negative, lead to inaccurate qualitative 
description of the atmospheric situation. Errors in the second component lead to imprecise 
computations of the optical and radiative parameters affected by layers. 

Errors in the determination of layer boundaries from lidar profiles are largely controlled by the 
signal to noise ratios of small signals. The crucial objective of the boundary algorithm is to find a 
threshold small enough to detect signals from rarefied layers but large enough to reject random 
noise as layers (or false positives). A large percentage of layers could be detected simply by using 
a single constant backscatter coefficient, based upon the computed molecular backscatter 
coefficient, as a threshold. If such a threshold were greater than the molecular value, the 
boundaries would be known to acceptable accuracy for many purposes. No false layers would 
result if the threshold were high enough. However, many significant rarefied, optically thin layers 
would be overlooked. If a threshold were too low, random noise would be often interpreted as layer 
signals. In both of these situations, the boundaries of identified layers would be uncertain. The 
occurrences of false negative results and false positive results are the competing detrimental 
effects in the selection of a proper threshold value. 
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The uncertainties associated with determination of layer boundary locations will be measured in 
terms of probabilities that boundary results are within specified confidence intervals. The 
magnitudes of these probabilities will be determined through studies of simulated profiles and 
proxy GLAS data. These studies will consist of application of the boundary algorithm to situations 
where the desired results are known. Comparison studies of the results of the output of the 
algorithm with the known situation will be conducted. Probabilities of deviations of the algorithm 
output from the truth will be computed from these studies and tabulated. Tables 6 and 7 illustrate 
two important types of confidence relationships that will be generated. 

Probability of layer detection failure 



fxil 



Ti 

Pi i 

P21 

P31 

t 2 

P 12 

P22 

P 3 2 

t 3 

P 13 

P23 

P33 


Table 6. Probability that the GLAS layer boundary algorithm will fail to detect an actual layer. T n 
represents threshold values, Hfcpresents the maximum backscatter coefficient in a layer and P mn 
represents probability of failure. 


Probability of layer boundary height error 
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Table 7. Probability that an analyzed boundary height will deviate from the actual boundary height. 
T n represents threshold values, Ah m represents a magnitude of height deviation and Pmn represents 
probability of failure 

Many types of relationships could be developed but those relating the layer/no layer result and the 
location of layer edges to selected threshold value are those that are most appropriate for the 
GLAS product output. 

3.3.4 Sample of Results 

The following figures are examples of results of the application of the layer boundary algorithm to 
simulated GLAS data derived from actual atmospheric conditions observed by the NASA ER-2 
Cloud Lidar System and to actual GLAS data. Figure 3.3.4 shows the analyzed layer heights of a 
single profile together with the threshhold profile and molecular profile. Figure 3.3.5 shows a typical 
example of cloud detection from actual GLAS data. Yellow hatch marks in the upper panel indicate 
where the boundary algorithm has determined cloud top while the purple hatch marks indicate true 
cloud bottom. Notice in Figure 3.3.5 that the true cloud bottom is only indicated for regions of the 
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Height AiVSl (km) Height (km) 


cloud where either a lower cloud layer can be identified or the ground return can be seen. Thus, 
much of the cloud between -56.1 and -42 latitude does not have a true cloud bottom indicated. The 
lower panel shows all cloud boundaries found, whether true or not. 


plot_cld_bsc3_10a.gif 



Figure 3.3.4. Simulated GLAS signal profile. Also shown are the molecular profile, the 
computed threshhold profile, and analyzed layer boundaries indicated with the symbols T 
and B. 
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Figure 3.3.5. Cloud boundaries obtained from actual GLAS data. In the upper panel 
cloud top is superimposed on the backscatter image with yellow dashes. True cloud 
bottom is denoted by purple dashes. Note that in the upper panel, where the signal is 
totally attenuated the true cloud bottom is u nkn own and is not plotted. 
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Figure 3.3.6. The output of the 1 064 cloud detection routine at the 4 second resolution. 
Cloud top is denoted by yellow dashes and cloud bottom by red dashes. 



Figure 3.3.7 An example of the 1064 nm 40 Hz cloud detection. The upper panel is the 1064 nm 
attenuated backscatter for a typical cloud scene and the bottom panel show the cloud layer tops 
obtained from an analysis of that data at 40 Hz. 
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3.4 Planetary Boundary Layer (GLA08) 


The height of the Planetary Boundary Layer (PBL) is one of the most important parameters retrieved 
from the GLAS atmospheric lidar data. PBL height is related to the fluxes of heat and moisture at the 
surface and can be used to estimate the bulk water vapor content of the PBL (Palm et al, 1 998). 
Because of the large aerosol gradient normally encountered at the top of the PBL, it is relatively easy 
to find the height of the PBL from high signal to noise (S/N) lidar data. The higher aerosol and 
moisture content of the PBL results in a much larger backscatter cross section, resulting in increased 
return signal. A strong inversion normally present at the PBL top traps the aerosol and moisture, 
thereby maintaining the large gradient of moisture and aerosol at the PBL top. The ability to measure 
the height of the PBL with both ground based and airborne lidar is well documented. Algorithms used 
with both types of data basically search the lidar signal for the large gradient of aerosol scattering 
within certain pre-defined levels of the atmosphere. Comparison of the lidar derived PBL heights with 
coincident radiosonde or dropsonde data has verified the accuracy of these methods. (Boers et al, 

1 984, 1 986; Melfi et al, 1 986; Palm et al, 1 998). 

Airborne lidars have frequently been used to gather high resolution measurements of tropospheric 
clouds and PBL structure over large areas (Melfi et al, 1986; Boers et al, 1991). Most airborne lidar 
systems consist of relatively large and powerful lasers which fly in the lower or mid troposphere. 
Consequently, the signal to noise ratio is high which makes the task of retrieving PBL and aerosol 
layer height from the lidar data fairly easy. The Cloud and Aerosol Lidar System (CALS), developed at 
NASA Goddard Space Flight Center, is an exception since it flies in the lower stratosphere and utilizes 
a relatively low power laser. Through the analysis of data from CALS and more recently, simulated 
GLAS data, we have developed schemes to retrieve PBL height from data with very low S/N (Palm 
and Spinhirne, 1987; 1998). This technique is described in section 3.4.11. 

Elevated aerosol layers (EAL) are not as ubiquitous as the planetary boundary layer, occurring only 
sporadically at various altitudes throughout the troposphere and lower stratosphere. Lidar is one of, if 
not the only remote sensing technique which can accurately resolve the height distribution of EALs. 
They are important because of their effect on the radiation balance and their contamination effect on 
many passive remote sensing measurements. The detection of EAL from lidar data is similar to that 
for the PBL height, but requires a somewhat different approach. Because of this, it will be addressed 
separately in section 3.4.12. 

3.4.1 Theoretical Description 

3. 4. 1.1 Planetary Boundary Layer 

Retrieving PBL height from the GLAS data can be difficult especially if the PBL is relatively dry and 
aerosol free. Even under the best of conditions (optically dense PBL and after sunset) it is very difficult 
to detect the PBL top from GLAS data on a shot to shot basis. Averaging of lidar shots to increase S/N 
is required. The degree of averaging depends on the optical depth of the PBL and lighting conditions 
(background noise). Linder typical conditions, the PBL top can be recovered after averaging between 
5 and 10 lidar returns. GLA08 is designed to detect the PBL height at two horizontal resolutions - high 
resolution (5 Hz or 8 shot average) which corresponds to 1 .4 km, and low resolution (1/4 Hz or 160 
shot average) which is about 30 km. There will undoubtedly be times when very little aerosol exists 
within the PBL, making the height determination very difficult or impossible at the high resolution. We 


51 



believe that at the lower horizontal resolution, we should be able to detect the PBL top well over 90 
percent of the time. There will also be times when ambiguities exist that tend to cloud the exact 
definition of PBL height (as defined by the lidar data). An example of this is when an elevated aerosol 
layer is riding directly on top of the PBL. In that case, it may be hard to discern the actual PBL top as 
distinct from the top of the elevated aerosol layer. 

GLA08 will use the 5 Hz, 532 nm attenuated backscatter profiles which are output from GLA07 for the 
calculation of PBL height. The algorithm must be designed to remove bad lidar shots and spurious 
noise spikes within shots. Failing to do so could result in noise spikes that are mistaken for PBL top. 
The filtering process can be done most efficiently by examining the quality flags that are output from 
GLA07. 

The PBL height algorithm processes the data in roughly 150 km chunks, which corresponds to 20 
seconds of data. The overall procedure is to first average 20 seconds of data to form one profile. That 
profile is searched below 7 km for the presence of the PBL and a ground return. If the PBL top is not 
found from this average profile, then it is assumed that the PBL top is not detectable for this segment 
of data and all the PBL heights for that time segment are set to zero. This would mean that the 100, 5 
Hz (high resolution) and the 5, % Hz (low resolution) PBL heights would all be set to zero. This is only 
expected to happen in cases where overlying clouds have attenuated the lidar beam, or in rare cases 
where the PBL is exceptionally devoid of aerosol. Now, there are certain criteria placed on the data 
within the 20 second data segment. First, if a cloud was detected for that shot (shot here means a 
single 5 Hz profile) via GLA09 above 5 km and the ground return was not detected, then that shot 
cannot be used in the 100 shot average. Further, if more than 50 percent of the shots fall into this 
category, then all the PBL heights for that segment are set to -1 . If a time gap of greater than 5 
seconds occurs, while forming the 20 second average, the 20 second average will have to be re- 
computed beginning after the time gap and all the PBL heights up to the time gap set to -2. 

Assuming that a 20 second average is successfully formed and that an average PBL height is 
detected, the next step is to go back through the 20 seconds of data and form five, 4 second (20 5 Hz 
profiles) averages and search each for the PBL top, using the 20 second average PBL top as a guide 
to where to search for the low resolution top. Similarly, when a PBL top is found from the 4 second 
average, the 20 5 Hz profiles that make up that segment will be examined individually for the high 
resolution PBL top, using as a guide the location of the 4 second PBL top. The output from this step 
represents the high resolution, 5 Hz PBL height. Thus, the general idea of the algorithm is to locate 
the PBL top at low horizontal resolution and gradually increase the resolution in a three step process. 
The exact technique used to locate the PBL top at any given resolution is discussed below. 

We need to identify the average ground bin (Gb) for the data segment under consideration. The 
position of the ground bin should not change substantially within a high resolution segment (5 Hz), but 
may change for a low resolution segment (4 seconds). For the 20 second average segment, the 
position of the ground bin could change substantially over mountainous terrain. The ground bin 
together with the last 20 second average PBL height in meters (H 20 ) gives us a reference from which 
to calculate various signal levels required by the algorithm. GLA09 will locate the ground bin from the 
532 nm return signal. When available, this will be used by GLA08 for the ground bin. However, there 
will be times when clouds attenuate the signal and no ground return is found. In this case, a calculated 
value of the ground bin will be used. Next, we need to compute the average signal level within the 
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boundary layer and above the PBL (within the troposphere). Let us call these average signals p pb i and 
Ptrop, respectively. We also need to find the maximum signal within the PBL. Let us denote this as 
pmax. The above filtering and averaging procedure should have eliminated all shots with no ground 
return and a cloud above 5 km. The reason that we do not want to eliminate all data with no ground 
return is that to do so would be to eliminate all cloud-capped boundary layer data. Instead, we want to 
eliminate all data with no ground return that was due to attenuation of the laser beam from mid and 
upper layer clouds, not from clouds that are associated with the PBL top. 

We begin by applying a 3 point binomial filter to the attenuated backscatter data below 7 km to form a 
smoothed profile (p s ): 


(3.4.1) p s (i) = S(l)P(i-l) + S(2)P(i) + S(3)P(i + l) for G b - 91 < I < G b 

where I represents the lidar bin number, G b - 91 represents the lidar bin corresponding to 7 km above 
the ground and S(j) is the binomial filter function with values: S(1 ) = 0.25, S(2) = 0.50, S(3) = 0.25. 

To obtain the average signal within the PBL (p pb i), compute the bin number that corresponds to half 
the average PBL height as k = H2o/(2.0*76.8). Then define the average PBL signal as: 


(3.4.2) p = (P (k - 1) + p (k) + p (k + 1)) / 3.0 

pbl s s s 

Similarly, to define the average signal above the PBL in the free troposphere (Ptrop), we compute the 
bin number that corresponds to 500 meters above the average PBL height as I = (H2o+500)/76.8. 
Where I is constrained to be greater than G b — 91 . The average signal above the PBL is then: 


(3.4.3) p = (p (l - 1) + p (/) + p (l + 1)) / 3.0 

trop s s s 

Next, define a signal level (Pt): 

(3.4.4) p t = fi trop + F phl (P pbl - p trop ) 

where F pb i is a threshold factor between 0.0 and 1 .0. In practice, the value of F pb i may vary from about 
0.4 to 0.7. A discussion of how to estimate the magnitude of F pb i is given in section 4.3.1 . Finally, we 
find the maximum signal between bins k and I. Call this p m ax, occurring at bin m. The algorithm then 
searches from that point (bin m) upward until 2 consecutive bins have signal values less than pt. The 
lidar bin corresponding to the top of the PBL is considered to be the first bin that is less than Pt. If we 
call this bin n, then the height in meters above ground of the PBL is: 

(3.4.5) H = (G - n) 76.8 

Pbl b 
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An example of a typical GLAS return for a clear marine boundary layer is shown in Figure 3.4.1 . The 
increase in signal due to the trapped moisture and aerosol within the boundary layer occurs at about 
900 m in this case. The various signal levels discussed above are labeled on the figure. 



0.0000 0.0010 0.0020 0.0030 0.0040 

Backscatter Cross Section (1 /km— sr) 

Figure 3.4.1. A nighttime simulated GLAS lidar return at 5 HZ showing the increase in signal 
associated with the marine boundary layer (below 1 km) and the various signal levels that would be 
computed by the algorithm from equations 3.4.2, 3.4.3 and 3.4.4 . The threshold value ptwas 
computed with F P bi = 0.5. 

Note that H 20 is the average PBL height defined by the processing of the last 20 seconds of data. This 
means that when we begin processing or resume processing after a large data gap, the initial value of 
H 20 must be assumed. While this is somewhat of a problem, it can be overcome by using the height of 
the maximum signal from the initial 20 second averaged profile as an estimate of H 20 . The maximum 
signal would be computed based on the data form 7 km altitude to 2 bins above the ground bin. 

After we have computed H 2 ofrom the 20 second average using the above procedure, we go back into 
that segment and form five, 4 second averages (20 shots). Each of these five profiles is searched for 
the PBL top in exactly the same way as described above, except for the following: the limits within 
which to search for the PBL top are more narrow. Now we use n - 5 and n + 5, which is a 750 m wide 
window centered on the 20 second average PBL height (H20). After each of these segments have 
been processed to obtain the low resolution PBL height (H4), the 20 shots which comprise them are 
individually searched for the PBL top in a similar manner, except we use a 600 m wide window 
centered on the low resolution height for that segment (H4). 
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Figure 3.4.2 Example of the output of the PBL retrieval algorithm using actual GLAS data. The red 
points in the bottom panel are the PBL top. Yellow denotes elevated aerosol. 
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Figure 3.4.3 Example of the cloud/aerosol discrimination routine. In the bottom panel, yellow 
denotes elevated aerosol while blue indicates cloud. 
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where F2 has the same value as in Equation 3.4.17. A maximum limit is also placed on p a as in 
equation 3.4.18, but now z corresponds to the height of the 20 second aerosol layer bottom. The 
search begins at bin b2+12 and continues to bin b2-12, searching upwards again looking for 2 
consecutive bins with a signal value greater than p a . The first of the two bins above the threshold level 
(p a ) defines the bottom of the layer. The above process is repeated for each of the five, 4 second 
average profiles yielding the high resolution elevated aerosol layer heights for the 20 second data 
segment. 



Figure 3.4.4. The GLAS derived PBL height for October, 2003 (upper panel) and the ECMWF PBL 
height for the same time period (lower panel). 

3.4.2 Error Quantification 

The accuracy of PBL or elevated aerosol layer height retrieval is governed by a number of factors. 
These are: 

1 . ) Signal to noise ratio of the data (determined by background condition and laser energy 

2. ) Sampling frequency (bandwidth, which determines the vertical resolution of the data). 

3. ) Number of lidar shots averaged together (horizontal resolution) 

4. ) Optical depth of the PBL 

Factor one encompasses numbers 4 and 5 as the signal to noise ratio increases with the square root 
of the number of shots averaged and the optical depth of the layer. Thus, the ability and overall 
accuracy with which we can detect the PBL top at low resolution (30 km, or 160 shot average) is going 
to be much better than high resolution (1 .4 km, or 8 shot average). Since the sampling frequency is 
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every 76.8 meters in the vertical, under the best of conditions, with high signal to noise levels, the PBL 
or aerosol layer can be resolved to a vertical accuracy of +/- 76.8 meters. As signal to noise values 
decrease (corresponding to either higher horizontal resolution or smaller aerosol backscatter) the 
retrieval accuracy will diminish. The magnitude of this effect must be determined by applying the 
algorithm to simulated GLAS data. For the horizontal resolution we hope to obtain with the GLAS 
measurements (1 .4 km), 76.8 meter vertical accuracy is too optimistic. We estimate that for the high 
horizontal resolution, the PBL top can be retrieved to within 150 meters. However, under typical 
conditions we estimate that we can obtain average PBL and aerosol layer height to within 76.8 meters 
at low resolution (30 km). 

PBL height is normally defined thermodynamically based on a rapid increase of potential temperature 
with height (temperature inversion) and a simultaneous decrease in relative humidity. A number of 
studies have shown that a rapid decrease in backscatter cross section also occurs at the inversion 
height, allowing lidar to provide a consistent and accurate measure of PBL height. Thus, we expect 
that the technique described here will usually yield the thermodynamic height of the PBL, but there are 
times when this might not be true. For example, elevated aerosol layers can lay directly on top of the 
boundary layer or there are cases where the boundary layer is growing into a residual boundary layer 
from the previous day. In these instances, the gradient of scattering at the PBL top might be relatively 
small and difficult to detect. Instead of yielding the PBL height, the algorithm might pick up the height 
of the elevated aerosol layer or residual boundary layer. Unfortunately, these types of errors are 
unavoidable when processing lidar data autonomously, without human interaction. 

GLA08 calculates the height of the PBL or elevated aerosol layer with respect to the ground return bin. 
If the ground return has not been detected from the data (from GLA07 or GLA09 processing), then 
clouds are assumed to be present. In this case, the algorithm will rely on either the last ground return 
bin found or a calculated ground return bin based on the onboard DEM value used for that shot. The 
onboard DEM is accurate to about 200 meters. Thus, when there are clouds obscuring the ground 
return, there could be a 200 meter error in the determination of aerosol or PBL layer top. 

It should be noted that the vertical resolution of the data combined with the technique used to find the 
PBL height place an upper and lower bound on the height of the PBL which can be detected by the 
PBL algorithm. It is estimated that the algorithm will have trouble detecting boundary layers which are 
less than about 150 m thick. Boundary layers this thin are relatively infrequent but do occur at times 
near the center of subtropical high pressure ridges over the oceans (and possibly elsewhere). 

Similarly, the algorithm as coded will not detect boundary layers that are higher than 7 km. As far as 
we know, the highest boundary layers occur over the Sahara desert and normally range from 4 to 6 
km in height. 

3.4.3 Confidence Flags 

Confidence flags for both the PBL height and elevated aerosol layer height can be constructed out of 
the difference between the average signal levels outside of the layer and inside the layer (ASI). These 
levels are computed by the respective algorithms as detailed in sections 3.4.1 .1 and 3.4.1 .2. Basically, 
the confidence in our height determination increases as ASI increases. If p aer and p a bv represent the 
average signal levels inside and above an aerosol layer (or boundary layer) respectively, then we can 
form the ratio: 
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(3.4.21) a = {P aer ~P abv )ip abv =X SI/p abv 

If the value of a is less than or equal to zero, then there is no confidence in the resulting height. As a 
increases, so does the confidence in the corresponding height measurement. It should suffice to 
compute and store a for the confidence levels of both the PBL height and elevated aerosol layer 
height. 

Another measure of confidence that could be used is the standard deviation of the heights for a given 
segment. Normally, for segments less than a few hundred kilometers, the PBL heights should have a 
standard deviation on the order of 200 to 400 meters. Any significant deviation outside of this range 
may indicate trouble with the algorithm. This approach could also be used for the elevated aerosols, 
except that the standard deviation is expected to be somewhat less, perhaps 50 to 200 meters. 

3.5 Blowing Snow Detection 

Blowing snow detection was not originally envisioned to be part of the GLAS atmospheric data 
products prior to launch. The idea for it came about after launch when inspection of backscatter 
images over Antarctica revealed what looked to be regions of thicker and brighter ground returns. 
Upon closer inspection, it was determined that these areas were thin (100 - 200 m) scattering 
layers in contact with the ground. Further investigation showed them also to be related to areas of 
high wind speed. Digging still further, the altimetry data in these areas often showed the 
characteristic range delay associated with multiple scattering which is most pronounced from low 
scattering layers such as blowing snow. Realizing its potential importance to the altimetry mission, 
blowing snow detection became a priority. Blowing snow parameters are stored on the GLA09 data 
product and include layer height and optical depth 

The blowing snow detection algorithm interrogates the lidar return bins directly above the ground 
for an elevated backscatter signal indicative of a scattering layer in contact with the ground. In 
order to accomplish this, it is imperative that the ground return bin be located so that it is certain 
that we are indeed looking at the bins immediately above the ground and that any contamination 
due to the ground signal itself be eliminated. The first step is to use the 1x1 degree Digital 
Elevation Model (DEM) to define a 400 m wide window within which a search for the ground return 
is performed. The ground return is generally the strongest signal in the lidar backscattered return 
provided that overlying particulate layers have not strongly attenuated the laser beam. The ground 
search is performed from 200 m below the (DEM indicated) ground working upwards in search of a 
signal exceeding the ground signal threshold (1.0x1 O' 3 nr 1 sr 1 ). The ground threshold was 
determined empirically and represents a value attained when the atmospheric column two-way 
transmission is greater than about 1 0%. If the ground signal is found (call this bin “G”), and If the 
1 0 meter wind speed at the current location is greater than 5 m s- 1 and the backscatter signal in the 
bin immediately above the ground (G-1) exceeds the blowing snow threshold (about 2.5x1 O' 5 nr 1 
sr 1 ), then a low-level, wind-induced “scattering layer” is assumed to be present. 

The blowing snow threshold is constructed from a scaling factor times the magnitude of 532 nm 
attenuated molecular (Rayleigh) scattering at the height of the current retrieval location. The 
scaling factor has a value of 20.0 and was determined by an iterative approach of adjustment and 
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review of retrieval results until they were satisfactory. The resulting threshold must be great 
enough to insure minimal false positive detections while small enough to retain adequate 
sensitivity. The algorithm then interrogates the bins above (G-2, G-3, etc.) until the backscatter 
within the bin is less than the blowing snow threshold backscatter level times 0.20 (about 5.0x1 O' 6 
nr 1 sr 1 ). The top of the scattering layer is then the last bin to exceed this threshold. Once the layer 
is defined, the gradient of backscatter within the layer is computed. If the gradient decreases 
upward, then the layer is assumed to be blowing snow. Conversely, if the backscatter gradient 
increases upward, the layer is assumed to be low fog or cloud. This check, designed to eliminate 
low clouds, is done only if the layer is thicker than 1 bin. The optical depth of the layer is then 
computed using an assumed extinction to backscatter ratio of 20 sr. 

The output of the blowing snow detection algorithm has been extensively checked for consistency 
and quality by generating and reviewing hundreds of images of the detected blowing snow layers. 

A limitation of the lidar technique is that the blowing snow layer has to be at least 50 m thick in 
order for enough backscatter signal to be collected in the bin immediately above the ground. This 
means that shallow blowing snow layers, which may be frequent, will probably not be detected. 
Further, blowing snow cannot be detected beneath thick or highly attenuating layers (tropospheric 
or polar stratospheric clouds with optical depth > about 1 .5), since detection of a strong ground 
return is required. The latter limitation implies that most of the blowing snow associated with winter 
storms (cyclones) will go undetected. These limitations will certainly result in lower blowing snow 
frequencies than actually exist. Furthermore, the magnitude of the discrepancy will depend on the 
cloud cover frequency of a given region. For instance, along the coast of Antarctica where blowing 
snow frequency is known to be high, it is also cloudier than more inland regions. 



Figure 3.5.1 Blowing snow frequency for October 2003 as derived from analysis of ICESat data 
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3.6 Optical Properties of Cloud and Aerosol Layers (GLA10 and 11) 

Before we examine the equations that will be used to retrieve optical properties of clouds and 
aerosols from the GLAS atmospheric data, we will present a short description of the physical 
processes which govern the light-particle interactions and the notable difficulties in using these. 

First, we note that the primary atmospheric observation channel of GLAS will be at 532 nm. Gas 
absorption processes are negligible compared to scattering processes at that wavelength and so 
they will be omitted in the derivation of the particulate optical depths given herein. Ozone 
absorption, although small, will be factored out prior to the optical properties algorithm by 
calculating ozone transmission profiles from ozone climatological databases then dividing the lidar 
signal profile (attenuated backscatter ) by the ozone transmission (See section 3.2). 

The observed or effective optical depth of a horizontal layer of particles between the orbiting lidar 
and a given altitude is the logarithm of the ratio of a laser’s initial normalized pulse energy to its 
energy at that altitude. Thus, the basic physical effect which permits finding the cloud and aerosol 
layers’ optical depths is the diminution of the lidar pulse energy as it is scattered or absorbed by 
the atmospheric constituents. As the laser pulse travels from the instrument, its photons interact 
with and are scattered (redirected) by the molecular and aerosol particles of the atmosphere. The 
lidar detector measures those photons which are redirected into a small solid angle centered at 
180 degrees (backscattered) into its receiver. The number of laser pulse photons received in a 
short time interval from a given atmospheric volume are recorded. This quantity is the lidar signal 
strength. It is proportional to the densities of particles in that volume and the combined scattering 
characteristics of the particles. These scattering characteristics are strongly affected by the shape 
of the scattering particles and the size of the particles relative to the wavelength of the laser light. 

A given scattering volume may contain zero (in a vacuum) to several scattering species each of 
which has its own density, size distribution, and scattering characteristics. 

A major challenge in optical analysis of lidar signals from cirrus clouds is that these clouds are 
composed of particles whose shapes and size distributions are not readily discernible by any 
remote sensing techniques. This forces us to incorporate some crucial assumptions in order to 
obtain quantified results. The validity of these strongly rely upon former experience with cirrus lidar 
observations (Spinhirne et al., 1990,1996). 

In particular, when attempting to obtain cloud optical depth from a spacecraft lidar, two 
assumptions are required regarding the scattering characteristics of the cloud. One of these 
assumptions is that the multiple scattering effect can be reliably quantified. The multiple scattering 
effect is the modification from the true optical depth caused by the increase in detected signal 
strength due to the portion of the detected signal which has experienced more than one scattering 
interaction. It is primarily the result of photons that are deflected only slightly during the scattering 
process. This is referred to as forward scattering and it serves to decrease the perceived optical 
depth. Ice particles typically have a very pronounced forward scattering component which will 
cause the multiple scattering effect to be quite significant. Multiple scatter is also a factor for 
aerosols, though much smaller. From calculations (Spinhirne, 1982), it is estimated that aerosol 
multiple scattered signals will have less than an 8 percent effect for even the most hazy conditions. 
Details of the procedure to handle multiple scattering are discussed in section 3.6. The other 
assumption is that the value of the extinction to backscatter ratio is known. The extinction to 
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backscatter ratio is the total scattered energy divided by the amount of backscattered energy. For a 
given scattering layer, it is assumed to be constant. Sometimes, under favorable circumstances, 
this ratio can be estimated from remotely sensed data, but computations based upon satellite 
observations often will require externally computed values. These are discussed in sections 

3.5. 1.1 and 4.5.1. The values of both of these parameters are determined by the details of the 
volumetric scattering phase function that quantifies the scattering effect as a function of scattering 
angle. 

3.6.1 Theoretical Descriptions 

3. 6. 1.1 Transmittance Solution to the Lidar Equation and C alculation of Backscatter Profiles 
(GLA10) 

The goal of the optical properties analysis of the GLAS lidar signal is to obtain particulate extinction 
cross section profiles (<j p ) and particulate layer optical depths (x p ), involving clouds, aerosols, and 
polar stratospheric clouds (PSC’s). The discussion given below essentially restates a derivation 
given many times in the literature. For example, see Spinhirne (1980,1996), Elouragini (1995), and 
Marenco (1997). The derivation of the multiple scattering factor (r|) will be handled in its own 
section (3.6). At this point one needs to note that transmittance, extinction, and optical depths 
obtained directly from the solution of the lidar equation are actually the apparent or effective values 
(Platt, 1979) without multiple scattering effects factored out and are denoted with the superscript 
prime. 

The working lidar equation for a spaceborne and nadir pointing lidar has been stated previously 
and can be rewritten in the following form: 

(3.6.1) P n (z)= p(z)T' 2 ( z). 

The left side of the equation is the calibrated normalized lidar signal of attenuated backscatter 
coefficient corrected for ozone attenuation. The total (particulate and molecular) volumetric 
backscatter coefficient at distance z is denoted by /?(z) and the two-way particulate and molecular 

transmission factor from z to the spacecraft altitude is given by T' 2 (z) , expressed as 
exp[-2(r m (z) + r),(z))] . Optical depth is represented by the symbol x, while the subscripts m 

and p designate molecular and particulate contributions, respectively. Furthermore, the influence 
of the multiple scattering effect (r|) on the particulate optical depth is described by: 

(3.6.2) r’ p (z) = jcr; ( z)dz =| r/(z)cj p ( z)dz = pj ( z)dz , 

where a p is particulate extinction. 

Since the molecular contribution to the total backscatter and transmission can be computed from 
theory, it is advantageous to separate the scattering terms into components which represent the 
molecular and particulate contributions independently. 

with p = and r 2 = r p x 
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the lidar equation becomes: 

( 3 . 6 . 3 ) p n = p p r P X+PJ' P X- 

The following relationships must now be defined: 

(3.6.4) T' p 2 =e 2 X dz 

(3.6.5) T,;=e J 

where S' and S m are the effective particulate and molecular extinction to backscatter ratios, 
respectively. T 2 (z) can be calculated accurately given the vertical temperature and pressure 
structure of the atmosphere from MET data or appropriate standard atmosphere data and the fact 
that Sm is known to be 8 tt/ 3 throughout the vertical profile. The purpose of this derivation is to solve 
the equation for the vertical profiles of fd p . The true particulate optical depth and extinction profiles 

can be then be computed from the values of S p , j3 p , and 


and S ' P = ~^~ ( assumed a constant for each layer), and 
and S m = , 


From these relationships, we see that: 


(3.6.6) 


dz 


= r 2 (-2S'k 


We can use this relationship to substitute for p p in (3.6.3) to arrive at: 

= 4k 2 ) 


(3.6.7) P„=-| 

4n 2 ) 


v/ 2S K 


dz 


■+ B T 2 T' 2 

r' m m p 


or 


dz 


2 s' P 

-2 S'BT' 2 =-^X 


By specifying z as the independent variable and T' 2 as the dependent variable, this is a first order 

linear ordinary differential equation; it is a special form of the Bernoulli equation. The solution can 
be found by using the common integrating factor method where the integrating factor is 

-2X f <j dz S . 

F = e J , and X = — . The general solution is: 

S L . 


(3.6.8) T 2X T' 2 = -2S’ p \ T 2{x - l) P n dz + K 

where the integrand is defined only where particulates are present and K is a constant of 
integration. 
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For convenience, we define the coordinate z as the vertical distance from the spacecraft, 
increasing downward. At this point we must now allow for the effect of the lidar pointing off-nadir at 
a zenith angle of 0. If we visualize the situation where the lidar pulse encounters layers of 
particulates after traveling through the molecular atmosphere from the spacecraft, we can define 
the boundary condition at the top of any layer, I B (z t ) , as: 

(3.6.9) I B (z t ) = T; 2sec0 ( Zi )T* Xsecd (z t ) , 

where Zf is the distance to the top of the layer. If the layer is the first layer encountered, the 
T' p 2sccd ( Zt )iem can be estimated as 1.00. The calculation of I B (z t ) for multiple layers is 
covered in detail in Section 4.5.2. 

So in general, the two-way particulate transmission within the particulate layer, whether cloud or 
aerosol, given a lidar zenith angle of 0 is 


(3.6.10) ^ . 

Ah \ Z ) 

This forward inversion processing continues throughout each particulate layer until T * ec9 p (z) < Tl or 
the signal from the earth’s surface is detected. Tl is a limit defined through error consideration 
(see section 3.6.2). Extensive automated use of this algorithm has been incorporated into the 
Global Backscatter Experiment (GLOBE) with aircraft lidar and into the Atmospheric Radiation 
Measurement (ARM) program with the ground-based Micro pulse lidar (MPL) with good results 
(Hlavka, 1998). Backward inversion processing, where the boundary conditions are known at the 
base of the layer, will optionally be used for low noise and high optical depth situations. Details of 
the backward inversion algorithm can be found in section 4.5.2 including equations 4.5.4, 4.5.5, 
4.5.6, and 4.5.7. 

An important ingredient of this transmission solution is the factor S p . When the particulate layer 
being analyzed is determined to meet the appropriate criteria for underlying signal analysis, an 
algorithm to calculate an estimate of S p will be called. If S pis found to be within tolerances, it will 
be used in equation (3.6.10). Appropriate criteria would be 1) layer is optically thin with either a 
lower layer or earth’s surface sensed and 2) enough clear air (no aerosols) exits below the layer to 
determine signal loss through the layer. The clear air zone must be at least a minimum thickness 
(around 1 km) and analysis is usually restricted to 3 km thickness. Ice clouds above 5 km are the 
most likely candidates. Under these conditions, an estimate of T p 2sec0 (z b ) (and thus an estimate 

of effective optical depth for the layer) can be found using the following equation, where zp is the 
distance to the bottom of the layer and z c is the distance to the end of the clear air analysis zone: 
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\P n (z)dz 

(3.6.11) r p 2sece (z b ) = g . 

T:r\z h i\p m {z)T;r ce {z)dz 

z b 

This method is called the integrated ratio technique. Simulations show this method is more stable 
under noisy conditions compared to other methods such as the log signal difference technique (see 
section 3.6.2). By defining I B {z b ) = T' 2sec9 (z h )T 2Xsec \z h ) , S' p can then be calculated through 
an iterative solution from the following equation: 

(3.6.12) S' p = 

2 sec 6\P n (z)T 2{X - l) ™\z)dz 

z t 

The iterative process is started with an initial guess of S' p as it relates to the X parameter, with the 

next iteration using the calculated value until the solution converges to a set tolerance. A version 
of this routine has worked well during automated MPL processing of aerosols using the calibrated 
signal to resolve the layer optical depth similar to the loss of signal in a cloud (Spinhirne, 1999). 

This routine should also function for PSCs and enhanced upper tropospheric aerosol layers. Later 
versions of this ATBD will look into the feasibility of expanding the integrated ratio technique by 
combining two close layers if they are the same layer type and the bottom layer meets the criteria. 
An average Sp can then be calculated to represent both layers. 

For atmospheric layers where S' cannot be calculated, a value will be assigned for each layer 

based on pre-defined look up matrices of S p and 77 , distinguishing between different cloud and 
aerosol regimes. Because the calculation of S' requires a clear air zone below the layer, all 

Planetary Boundary Layers (PBL) will have to default to the pre-defined matrices. Details of the 
decision matrices for S p look up tables are presented in section 4.5.1 . S p will be determined as: 

(3.6.13) S' p =T]Sp, 

where the multiple scattering factor, is separately estimated from appropriate look up selection 
distinguishing between apparent cloud or aerosol type, layer top and bottom heights, effective 
optical depth estimate, and particle size (see section 3.6). Initial determination of S p for clouds will 
be driven by cloud temperature. The underlying surface signal attenuation is an additional factor to 
improve retrievals. 
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Normalized Backscatter Phase Function Frequency of Occurrence 



Figure 3.6.1 Phase function (1/S P ) for Midlatitude Cirrus Observations 

A lidar study of mid-latitude cirrus (Eloranta, 1 999) indicates that although S p can vary by 30 or 
more, by far the highest frequency of occurrence is near 24 sr (refer to figure 3.6.1). Water clouds 
have a much lower variation. Determination of S p for polar stratospheric clouds will be handled as 
a special subset of the aerosol look up table because they have more of an aerosol origin than a 
water origin and will be processed at the aerosol time resolution. Determination of S p for regular 
aerosol will be driven by geographic location, layer height, relative humidity, and possibly surface 
signal attenuation analysis and wavelength ratios of solar reflectance at 532 and 1064 nm, with 
geographic location the most important factor. Geographic location can be channeled into three 
main aerosol regimes: continental, desert, and maritime (Ackermann, 1998) with functions relating 
the influence of relative humidity. Analysis of the GLOBE data set of 1 990 suggests that, on 
average, aerosol S p equals 28+5 srfor all height levels, even though there were distinct boundary 
layer and upper tropospheric layers with different sources. An example is shown in figure 3.6.2. 

Note that if T ipi (z) = 1 , which means that molecular scattering is negligible at all processing levels, 
the transmission equation for nadir pointing lidar reduces to: 


(3.6.14) T' p 2 (z) = \-2S' p \P„dz ' , 

z t 


which many times is sufficient for cirrus cloud analysis using a 1 064 nm channel. 
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Figure 3.6.2 Despite complicated vertical structure, the GLOBE project showed that S p did not 
vary appreciably in the vertical. 

Finally, in order to obtain the relative density for aerosol and cloud scattering, it is useful to solve 
the lidar return signal for the actual particulate backscatter cross section without attenuation. To 
solve for this backscatter cross section profile, use results from (3.6.10) as input to (3.6.3) by 
rearranging: 


(3.6.15) 


PM = ' 


PM 


C ( Z )V ( z ) 


P,M- 


3. 6. 1.2 Aerosol Extinction Cross Section 

Once the particulate effective transmission and backscatter profiles for each aerosol layer sensed 
have been calculated, it is a straightforward procedure to determine the aerosol extinction cross 
section profiles. Extinction cross section for particulates (a p ) is defined as the total scattered 
energy at height z or the change in optical depth with height: 

dr n (z) 

(3.6.16) a p (z)= p , where 
dz 
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r p is the particulate optical depth at distance z and is defined in Section 3.6.1 .4. However, since 
the backscatter cross section profile has, at this point, already been calculated, the aerosol 
extinction to backscatter ratio (S p ), earlier obtained through table look up or its relationship with S' p , 
could also be easily used to obtain the extinction profile for visible wavelengths: 

(3.6.17) o- p (z) = S p p p (z) 

Note that multiple scattering has already been accounted for in the calculation of S p . 


3. 6. 1.3 Cloud Extinction Cross Section 

As discussed in section 3.6.1 .1 , the calculation of the extinction cross section of clouds from 
backscatter lidar data requires knowledge of the 180 degree scattering phase function, or 
extinction to backscatter cross section, (S p ) and a correction, in the case of space born lidar 
especially, for multiple scattering (rj). In all cases an extinction solution, or correction, for cloud 
lidar can only be applied for a limited range of optical thickness. Our experience with ER-2 remote 
sensing indicates an upper limit of approximately 1 .5 effective optical depth. Signal to noise issues 
and others will be a factor for other systems, and modeling and testing with simulated GLAS data 
will determine the applicable limit. In order to determine the effective attenuation, neglecting first 
multiple scattering, most generally previous work has made the assumption of a constant phase 
function within cloud layers. With this assumption it is well know that integration of the observed 
attenuated backscatter cross section for optically thick clouds is equal to half the backscatter to 
extinction cross section 

(3.6.18) \B\z)dz=> kllrj as * => 00 

By identifying the limiting integral value, a solution for the effective backscatter to extinction value is 
known. For cirrus, Platt (1979) and others have used infrared emittance to determine asymptotic 
values. For nadir observations, Spinhirne and Hart (1 990) have shown that the disappearance of 
the surface signal below the cloud can be used to identify the asymptotic value. For real time 
processing however there are limitations. The assumption must be made that near by thin clouds 
are in character with dense clouds. Also signal noise and the complexity of real cloud formations 
can be expected to introduce significant error, based on ER-2 experience. The noise effects and 
appropriate application routines can be examined from modeling. A more basic limitation is that 
the multiple scattering correction to the derived effective extinction is already a large uncertainty 
term, and complex algorithm development for the effective attenuation may not be warranted. 

Another approach to obtaining the extinction cross sections, which is the one we prefer for 
automated processing, is to start with assumed 180 degree scattering phase functions. For water 
cloud this is accepted as a good assumption where 17.8 (sr) is widely applicable (Spinhirne et al., 

1 989; Pinnick et al., 1 983). For cirrus, modeling has not shown such an universal value, give the 
complexity and variation of cirrus shape and size. However experimental measurements have 
shown, likely because most cirrus are complexes of many different crystal types, that cirrus phase 
functions values tend to peak statistically toward characteristic values (E. Eloranta, personal 
communication; Spinhirne et al.,1996). With further work it will be possible to tailor values to 
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geographic location and cloud temperature or height. When the profile is determined to be 
appropriate for direct S ' p analysis, an algorithm to calculate an estimate of S ' p will be called and the 
calculated value will be compared to the table. Similarly, modeling for the multiple scattering 
correction, discussed in a later section, will lead to look up tables for the correction based on the 
cloud location and temperature. 

Given a 180 degree phase function value, we use equation (3.6.17) to get the vertical profile of 
visible extinction cross section, a p of a cirrus cloud from top to the upper limit of the effectiveness 
of an attenuation correction. Please note that multiple scattering has already been accounted for in 
the calculation of J3 P and S p . 

Any conversion of the visible extinction coefficients in clouds to thermal infrared absorption 
coefficients will rely on the assumption that the ratio of these two parameters is constant through 
the vertical extent of a cloud layer. A profile of lidar derived backscatter coefficients can be 
converted to absorption coefficients by a direct multiplication. The value of this constant absorption 
ratio, q, can be approximated from the results of theoretical studies. The investigation into cloud 
infrared absorption conversion will be investigated as a level 3 product. 

3. 6.1. 4 Cloud and Aerosol Layer Optical Depth (GLA11) 

A fundamental use of the spaceborne lidar is to detect and quantify the layers of optically dilute 
clouds that often reside in the mid to high troposphere where the temperatures are cold . These 
temperatures result in low water vapor density. Because the total amount of water must be 
conserved, when clouds form, the particle density will likewise be low. Clouds which form at lower 
altitudes are generally denser because of the greater availability of water. In such clouds, the 
useful geometrical penetration of the lidar signal is limited because of laser pulse attenuation. 

The clouds of interest are generally classified as cirrus clouds. They are usually composed 
primarily of ice crystals with some coexistent supercooled water droplets. Analysis of PSC’s are 
also of strong interest. Both these clouds often have a sufficiently small optical depth that a 
meaningful lidar signal can be detected at the bottom of the cloud layer. In these cases, a total 
layer optical depth can be derived by the lidar. Sometimes, two or more layers exist and the optical 
depth of each layer can be determined. 

The lidar signal can also be used to estimate the optical depth of the layer of non-cloud aerosols 
which reside in the planetary boundary layer. These aerosols can be composed of a variety of 
substances that are trapped by the temperature inversion which tends to exist at the top of the 
boundary layer. In some cases elevated haze layers of significant density are also found higher in 
the troposphere or stratosphere which have appreciable optical thickness. Examples include 
volcanic dust layers, smoke layers, and dust layers caused by periodic continental dust storms. 

The boundaries of any of these layers that are significant would be located by using the 
backscatter discontinuity algorithm of Section 3.4. 

The solution to the lidar equation to obtain particulate effective optical depths ( t’ ) at any range z 

from a nadir pointing spaceborne lidar is by definition a straightforward relationship with the 
particulate effective transmission calculated in (3.6.10): 
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(3.6.19) 


( i A 


;(z)=- - in[r; 2 (z)] or r p \z) = e 


-2r' (z) 




This solution will achieve best results in terms of the magnitude of error when applied to situations 
where the optical depth is relatively small. To evaluate and quantify this declaration we examine 
the relationships from which x' p is computed from cirrus data: 

We neglect molecular scattering within the cloud such that 

z 

(3.6.20) T' p 2 {z) = 1 - S' p 2y(z) where y(z) = J P n dz 

z t 

We see that T' 2 {z) approaches zero as 2y approaches 1/5" . Random noise excursions 

superimposed upon the detected signal can cause the computed value of T' 2 {z) to become less 

than 0 as the integral to evaluate gamma is numerically computed from the lidar signal. In this 
situation, x' p becomes undefined. Also, differentiation of (3.6.19) and (3.6.20) shows 

S'dy 

(3.6.21) dv' p =^f 


This means that a given excursion dy causes an error in r' p in inverse proportion to the value of 
T' 2 ; that is, the magnitude of the error becomes larger as the effective transmittance become 

small and the effective optical thickness becomes relatively large. Based upon experience gained 
from aircraft lidar studies (Spinhirne, 1990), computational errors in cloud optical depth for GLAS 
due to random noise remain tolerable until the value of T' 2 reaches 0.12-0.20 or r' p =1 .1-0.8. 

Where the clouds are more optically thick, the lidar cannot give meaningful results. We will discuss 
the details of computational uncertainty more fully in section 3.6.2. 

The specific method we will be using to calculate the particulate layer optical depth stems from the 
same transmission solution to the lidar equation, put uses the relationship of the extinction cross 
section profile in the layer (described in sections 3.6.1 .2 and 3.6.1 .3) to optical depth. The final 
optical depth products from these calculations will be the optical depth ( r, ) for each of the 
particulate layers meeting the analysis criteria: 

z b 

(3.6.22) t, = J <j p (z)dz , 

z < 


where z b and z t are the bottom and top locations of the particulate layer, respectively and 
multiple scattering has already been factored out. 

The vertical coordinate limits on the integration in the transmittance equation in (3.6.10) will be 
determined by the cloud and aerosol boundary algorithms described in Sections 3.3 and 3.4. In 
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practice, the integration will be carried out starting at the first particulate layer top. Although the 
whole molecular transmission vertical profile starting at the top of the atmosphere and ending at 
the bottom of the lowest particulate layer sensed will have to be calculated, the particulate 
transmission vertical profile will be calculated only inside cloud and aerosol layers. The boundary 
condition (equation 3.6.9) at the top of any secondary layer will involve the particulate transmission 
squared at the bottom of the layer above and the molecular transmission squared at the top of the 
current secondary layer. 

The attenuation of the pulse energy due to molecular scattering in the intervening clear air layers is 
small in the mid to high troposphere where the optically thin clouds reside. The magnitude of the 
molecular scattering is a significant fraction of the aerosol scattering since the gaseous 
atmosphere is relatively dense at the low altitudes of the boundary layer and the optical density of 
the aerosol particles are typically much lower than that found in even cirrus clouds. 

Optical parameters would be obtained either empirically or from prior studies for aerosol layers. In 
practice, three necessary conditions for determination of the boundary layer or elevated haze 
optical depth will be that: 1) the top of the layer is detected, 2) there is no obscuring cloud layer 
present, and 3) the earth’s surface or a lower particulate layer has been found by the lidar. 
Evidently, the most prominent source of uncertainty will be how closely the actual aerosols which 
are being observed match the characteristics of the assumed aerosol type. Re-processing from 
level 3 extinction-to-backscatter investigations will help reduce these uncertainties. 


3.6.2 Error Quantification 

The most important optical measurements derived from the lidar measured backscatter profiles are 
the total effective optical thickness and transmission of particulate layers which are fully penetrated 
by the laser pulse. We will inspect here the effect that various uncertainties have on the uncertainty 
of the derived values of these physical quantities. 

To begin, we will use the relationships 

(3.6.23) T' 2 (z b ) = 1 - 2 S' p y{z h ) (ignoring the molecular component), 

(3.6.24) t' p (z b ) = - ^ In \r' p 2 (z b )] , and define the parameter 

(3.6.25) a = S' p y(z b ) , where 0 < a < 0.5 and y = J P n dz 


The subscript p denotes particulate and a will represent actual value of the product which the 
measurements are attempting to attain. Because practical computations become unstable for 
relatively optically thick clouds, the useful limits of lidar measurements are exceeded as the value 
of a goes above 0.475. 

Differentiation gives us, 
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(3.6.26) dT ' p 2 = -2 ydS' p - 2 S' p dy = -2 


a- 


ds 

51 


- + da 


and 


dz P 2 t' 2 


dT 


f 2 


If we let dS' p and da represent deviations from the correct values of the respective parameters, 

then we can assess the effects that such deviations will have on the derived values of these 
parameters. To simplify this assessment, we will estimate the effects of each deviation 
independently. 



a 

Figure 3.6.3 Computed sensitivity in optical depth from error in S p ’ 

Figure 3.6.3 shows an example of expected error if da = 0. The plot in the figure summarizes the 
percent change produced in the computed values of z' p by an error of 25 percent in the estimate of 

S' p . The importance of such errors is determined by what the purpose of the computed values 
are. 

The plot in Figure 3.6.4 illustrates expected magnitude of deviation in the computed optical depth 
as a function of a when dS' p = 0 and there is a typical 5 percent error in the magnitude of the 

integrated backscatter. We see that the magnitude of the error in the optical thickness becomes 
very large as the limit in meaningful measurements is approached at a « 0.475. For larger errors 
in the evaluated magnitude of a , the uncertainty in z' p is even larger. A fact that reduces the 

detrimental effect revealed by these relationships is that, for a given evaluation of optical depth 
from a lidar profile, the random fluctuation contribution to dy will become smaller as more 
samples are used to compute the result. This means that for layers of a given optical depth, the 
error in the optical thickness will be less for layers of greater geometrical depth. These are typically 
the types of layers of cirrus and aerosols which are the greatest interest to climatological studies. 
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The accuracy of the retrieved backscatter lidar signals relies heavily on the signal to noise ratio of 
the data. The signal to noise ratio rises and falls with the following: 

1) inverse of the strength of the background signal, 

2) strength of cloud or aerosol return, and 

3) horizontal smoothing of lidar shots. 

Tests using simulated backscatter data developed from GLAS instrument specifications as of 
January 1 , 2000 were run to estimate the accuracy of lidar signal techniques for extracting the 
extinction to backscatter ratio directly from the signal return of elevated layers and to estimate the 
accuracy of optical depth retrievals as per the operational algorithm. These tests were performed 
using three different background lighting conditions: 1 ) no background, 2) daytime over dark ocean, 
and 3) daytime over bright cloud. Figure 3.6.5 shows simulation results comparing a layer with no 
background noise with the same layer after bright daytime noise was applied. Resultant retrievals 
using 1 second averaging (simulating cloud processing) show the drop in accuracy in S p and 
x calculations with increased background noise. For these simulations, r\ was set to 1 .0. 
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The effect of noise on GLAS optical properties retrievals 


SIMULATED GLAS BACKSCATTER PROFILE SIMULATED GLAS BACKSCATTER PROFILE 

NO NOISE CASE . 1 SEC DAY. BRIGHT CLOUD CASE . 1 SEC 




Figure 3.6.5 Simulations showing the effect of noise on GLAS optical properties retrievals from a 
single layer of optical depth 0.45 

Figure 3.6.6 compares the log signal difference technique and the integrated ratio technique for 
overall accuracy in calculating Spfrom the lidar backscatter signal inside elevated layers. In 
general, although the log signal difference technique is slightly more accurate during no noise 
situations, the integrated ratio technique, our current algorithm, is much more stable during noisy 
signals with acceptable error. Errors for both 4 second (aerosol) and 1 second (cloud) are shown. 
Figure 3.6.7 displays simulated error results for optical depth retrievals also as a function of noise 
using the current protocode algorithm. The simulations show that on average for single layer 
conditions of moderate optical depth, output error will be near 4 to 5 percent for low noise 
situations but rise during noisy conditions to 20 percent for clouds. Aerosol errors remain much 
more stable with increasing noise. These simulations can be thought of as a “best case” scenario 
since multiple layers, low optical depth situations, high optical depth situations, and the addition of 
multiple scattering all tend to increase the error. 
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Figure 3.6.6 Accuracy of lidar techniques for extracting extinction to backscatter ratio (S p ) from the 
signal profile as a function of noise using GLAS simulated backscatter profiles for a single layer of 
optical depth 0.45 

Errors in the transmittance solution due to lidar signal degradation and atmospheric molecular 
misrepresentation are discussed further in Section 3.2.2. 

Real time error analysis will accompany the optical properties processing. Representative error 
profiles of the signal for the various time resolutions will be developed from the standard deviation 
profiles superimposed on the original signal. Optical processing of the error profiles will allow for 
the calculation of error bars for calculated S p , retrieved optical depth, and backscatter and 
extinction profiles. Details of the methodology can be found in sections 4.5.1, 4.5.2, and 4.5.4. 

The graphs in this section represent an initial analysis into quantifying errors in optical products 
from the GLAS atmosphere channel. Error analysis is on-going and will result in more detailed 
projections with further protocode testing using simulated lidar returns. Because the particulate 
transmission is restricted to a value of Ti and above (.35 <h< .45) to keep the integration stable, 
effective accumulative particulate optical depth is restricted to « 0.9 or below, and true optical 
depth is restricted to roughly 3.0, depending on the value of r\. With cloud profiles averaged to 1 
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second or 7.5 km and aerosol profiles averaged to 4 seconds or 30 km, we believe optical depths 
can be calculated to within an error of 50% in the troposphere. 
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Figure 3.6.7 Accuracy of lidar inversion retrievals for layer optical depth as a function of noise 
using GLAS simulated backscatter profiles 


3.6.3 Confidence Flags 

Confidence flags for GLA1 0 and GLA1 1 will include a measure of quality for the following 
parameters per layer: aerosol backscatter cross section, cloud backscatter cross section, aerosol 
extinction cross section, cloud extinction cross section, cloud optical depth, elevated aerosol optical 
depth, and boundary layer optical depth. Polar stratospheric clouds are incorporated into the 
aerosol analysis. The quality rankings, ranging from 0 to 15, will come directly from the magnitude 
of error bars created from the real time error analysis described in Section 3.6.2. Use flags will 
accompany the quality flags. The optical depth and extinction use flags are designed to contain 
atmospheric layer type information, such as water versus ice cloud, or maritime versus continental 
aerosol, or tropospheric versus stratospheric aerosol. This information will be passed to the use 
flags from the S ratio default decision matrices. The backscatter use flag will contain information 
on signal saturation conditions per layer. See Section 4.5.4 for a further discussion on quality 
control. 
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3.7 Multiple Scattering Correction 

The multiple scattering factor is a complex function of particle scattering phase function and the 
vertical distribution of the scattering plus the field of view and the height of the lidar receiver. A 
procedure for the correction of the GLAS lidar signal for multiple scattering from cirrus and other 
optically thin clouds is presented in this section. Two methods have been developed for 
calculations of GLAS multiple scattering using cloud and aerosol models. One is a precise Monte 
Carlo radiative transfer model, and the second is a computational fast analytic approximation. 
However a precise radiative transfer calculation to account for the effects of the multiple scattered 
contribution is not practical for real time analysis, and approximation by semi-empirical methods is 
necessary. For initial level 2 processing, the value of the multiple scattering factor may be decision 
matrixed into a look up table based on parameterized calculations. The best method to ultimately 
correct for multiple scattering is a subject in development but a preliminary acceptable procedure 
can be described. 

3.7.1 Theoretical description 

As was indicated in the optical properties introduction, some of the scattered energy in an 
atmospheric layer (particularly clouds) will undergo additional scattering and reenter the lidar 
detector’s field of view. This multi-scattered energy is indistinguishable from once scattered energy. 
The multiple scattering produces an ambiguity in the interpretation of the lidar signal. The laser 
pulse at any level z appears to be more energetic than would be calculated from the simple optical 
thickness attenuation. 

The component of the lidar signal from multiple scattering arises chiefly due to the strong forward 
scatter component, or diffraction peak, where the forward scattered photons stay within the 
receiver FOV. The width of the forward scattering is strongly dependent on particle size and may 
be approximated by 

d 2 s *X 2 /n 2 r 2 

where 0 S is the width of the forward diffraction peak containing one-half of the scattered photons 
and r represents the particle radius. If none of the forward scattered photons were lost from the 
FOV of the receiver, higher order scattering is neglected, the transmission term in the lidar 
equation can be written as e- 2 ( T - 1/2T ) . This is equivalent to e' 2 ’ 1T where r|=0.5. In actuality neither 
assumptions just mentioned will necessarily hold. Only for the cases where the diameter of the 
receiver FOV foot print at the top of a scattering layer is sufficiently large, and width of the forward 
peak sufficiently narrow will the approximation be close. However for space borne lidar the FOV 
footprint is large compared to terrestrial systems, approximately 100 m diameter, and for cloud 
scattering, especially cirrus, the narrow forward peak is expected. Under these situations prior 
work with ground based and airborne lidar indicate that the constant r| value near 0.5 is an 
acceptable approximation. A constant value of r\ can be assumed and variation in the value from 
0.5 determined for a first order correction of higher order scattering and loss of photons from the 
receiver FOV. For aerosol, or extended ranges due to multiple cloud layers, a constant r| factor or 
values near 0.5 will not hold. In any case, values of r\ as a function of layer heights and 
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propagation depth, and parameterized cloud and aerosol particle models, can be calculated as a 
basis for look up tables for real time processing. 

To account for the multiple scattering affect, we assert then the transmittance and optical depths 
obtained from the solution to the lidar equation are apparent or effective values. For lidar cirrus 
studies Platt (1981) proposed an extension of the single-scattering lidar equation to account for 
multiple scattering by introducing the parameter r| to the equation: 

(3.7.1) r' p (z) = jc7 p (z')rj(z')dz' = j S p /3 p (z')j(z')dz' and 

z t z t 

(3.7.2) T' p =e T ' p 

The parameter ?i is the multiple scattering correction factor where 0 < r\ < 1 . The superscript 
prime is used to denote the effective value. Modeling studies have indicated that usefully accurate 
results can be obtained if a constant value of jj is used within the integrations for typical cirrus 
layers. We use this to obtain: 

(3.7.3) z' p {z)=S' p \p p {z'}lz' 


where S' p = S p ij is denoted the effective extinction to backscatter ratio. 


The multiple scattering factor is accurately calculated by Monte Carlo methods or approximated by 
analytic methods (Duda et al., 1999). As mention above, the r\ coefficient as a constant value is 
inaccurate to apply toward some aerosol or cloud layers. However, for the case of cirrus clouds (or 
other clouds) where the cloud particle sizes are much larger than the wavelength of the lidar, r\ is 
shown to be a property of the forward diffraction phase function and can be computed analytically. 
The question to be answered by a parameterization is the appropriate r\ factor. 

Starkov and Flesia (1998) have developed an analytic formula to compute the multiple scattering 
correction factor (rfi for optically thin cirrus clouds. They assumed that the atmosphere was 
divided into N arbitrary layers (layers 1 through N) with at least one clear-sky layer (layer 0) above 
cloud top. Letting the clear-sky atmospheric phase function equal po(0,R), the cloud phase 
function in layer i equal pi(0), and 0i equal the width of the forward diffraction peak for particles in 
layer i, the multiple scattering correction factor for the nth layer at range R can be computed from: 

n - 1 

(3.7.4) r/(R)r(R) = 

n = 0 

(3.7.5) = + Tn(R), 

n = 0 
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where ii(R) is the optical thickness of layer i at range R inside layer Li 

(3.7.6) t{r)= f ( <Jm(x) + <Ja(x))dx + [ <Jc(x)dx, 

JRi JRi 

(a m (R) and a a (R) are the molecular and aerosol extinction coefficients in the atmosphere, and 
dc(R) is the cloud extinction coefficient) 

and 

(3.7.7) 7/0 = 1 — 2k f" po(0,R) s i n Q dO, 

Jo p„(k) 

(3.7.8) r/i = 1 - 2k f p(0) Pn ^~^ sin 6 dO. 

Jo p„{k) 

Wiscombe (1977) notes that the backscattering amplitude pi(7t) (and similarly the ratio 
Pi(7t-e)/pi(7i)) is one of the most difficult Mie quantities to calculate accurately, and can vary over 
orders of magnitude for small changes in particle size. If spherical particles are used to compute 
Pi(7t), the backscattering amplitude can be integrated over a broad size distribution to make it a 
smoother function of particle size. Mishchenko et al (1997) have calculated the scattering of light 
from polydispersions of thin disks and oblate spheroids. Like ice spheres, size averaging in the 
distribution will smooth out the scattering phase function. If the particles are horizontally oriented, 
the phase function will have a strong peak at the 180° direction. Macke (1993) has shown using 
ray-tracing calculations that the backscattering phase function is dependent on the shape of ice 
crystals as well, ranging from highly peaked functions for crystals having parallel or perpendicular 
planes (columns) to approaching zero for hollow bullets. Recent calculations and observations, 
however, suggest that for most cirrus clouds, the ratio pi(7t-0)/pi(7t) from equations 3.7.7 and 3.7.8 
may be less variable than might be expected from Mie scattering theory. Nicolas et al (1997) have 
shown that for clouds with optical depths of one or larger, it may be possible to compute an 
effective backscattering coefficient that is an average of the scattering properties around 180°. 
Also, analysis of extinction to backscatter ratios in cirrus from high spectral resolution lidar data 
(Eloranta, 1999, personal communication) shows that the observed scattering from backscattering 
angles does not vary as much as theoretical calculations of pure ice crystal shapes. The relative 
invariance of the observed backscattering coefficients is mostly likely due to averaging effects from 
the different particle shapes and sizes found in cirrus. Therefore, for clouds with optical depths 
greater than unity, equations 3.7.7 and 3.7.8 can be approximated as 

(3.7.9) 770 = 1 — 2/rf poi^O ,R)peff{7r^)svn6 dO , 

JO 

(3.7.10) 77 / = 1 — 2/r f pi(&)peff(K)sin6> dO, 

J 0 
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where P e ff(7t) is the effective backscattering coefficient. Note that in equations 3.7.9 and 3.7.10, if 
Peff(7t) is equal to one, then l-iq equals the portion of energy scattered in the forward diffraction 
peak. 

From diffraction theory, the width of the diffraction peak may be alternatively defined as 1 .21 XI d, 
where X is the lidar wavelength and d is the particle diameter. Using this definition of the diffraction 
peak width, the portion of the energy scattered in the peak can be calculated for ice spheres from 
Mie theory. The results are presented for the 0.532 nm channel in Figure 3.7.1 as a function of 
particle radius. For monodisperse spheres, the scattering in the diffraction peak oscillates. The 
central line in Figure 3.7.1 represents the diffraction peak scattering for a broad size distribution of 
particles, in which the size averaging tends to smooth out the oscillations. 


Diffraction peak scattering probability 



Figure 3.7.1 The portion of the energy scattered in the diffraction peak as a function of particle 
radius for ice spheres. For large particles, the portion approaches 0.42. 


As the particle size increases, the fraction of the energy scattered into the diffraction peak 
approaches 0.42. Nicolas et al. introduce a similar model where the amount of energy scattered 
within the forward peak (1 - iq) is given as 

(3.7.11) 1-77 = — + R, 
co 
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where co is the single-scattering albedo and R is a measure for nonspherical particles that 
describes the fraction of light refracted into the forward direction through opposite parallel faces. 

From the methods and results as described above plus other available knowledge, appropriate 
values to apply in calculations can be obtained. Also values of rj for cirrus analysis can be 
parameterized based on the height of the cirrus layers and observed depth. As an example of the 
effectiveness of approximations, values of r) beneath cirrus are shown as determined from 
accurate Monte Carlo calculations in Figure 3.7.2. For a given depth below the cloud the r\ value 
is seen to be independent of optical depth as required. In addition the value immediately below the 
cloud has an r\ of approximately 0.4. The increase for depths below cloud base more than 2 km 
greater are not dramatic. 



Figure 3.7.2 Cloud height plays an important role in determining multiple scattering from ice clouds 
referenced from the ground. The results of Monte Carlo calculations of the apparent optical depth at 
the surface as a function of the true optical depth for differing cloud layers is shown. The lower the 
cloud, the lower the ratio of effective to true optical depth, or r\. 

A cloud classification based on cloud temperature, geographic location, cloud height, and 
integrated backscatter in the layer will be used to parameterize a systematic cloud multiple 
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scattering factor look up table. The table will be generated by systematic Monte Carlo calculations 
supplemented by analytic models. 

For aerosol, there is typically not the sharp forward scattering peak as there are for clouds and 
their larger particles. The approximation for a constant r\ with depth is not expected to hold as well 
as for clouds. However the initial Monte Carlo calculations of GLAS parameters for the aerosol 
multiple scattering indicate that the multiple scattered component of the lidar signal is no more than 
20%. Also most generally the more optically thick aerosol are concentrated in thin layers at the 
surface. The approximation of a constant r| for aerosol will be tested. It is expected that the errors 
will not be a dominant uncertainty for optical depth retrievals. A scene classification based on 
geographic location, integrated backscatter in the layer, and aerosol height distribution plus a 
systematic aerosol multiple scattering calculation look up table will be used for an r| factor. Errors 
in ii extracted from look up tables based on modeled results have yet to be formulated. 


3.7.2 The Multiple Scattering Algorithm 

The multiple scattering factor (r|) is a complex function of particle scattering phase function, the 
vertical distribution of the scattering, plus the field of view and the height of the lidar receiver. 
Because of forward scattering, the effect of multiple scattering is to have more of the lidar pulse 
making it through the cloud or aerosol layer than would occur by single scattering. If more of the 
signal makes it through the layer, then the calculated (or effective) optical depth is perceived to be 
smaller. Applying a multiple scattering factor to the effective optical results corrects the results 
back to single scattering conditions: 

(3.7.12) t = , where t is the effective optical depth. 

A procedure for the correction of the GLAS lidar signal for multiple scattering from cirrus and other 
optically thin clouds, plus aerosol layers is presented in this section. An additional procedure to 
calculate the range delay or pulse spreading of the lidar signal caused by longer trajectories of 
multiple-scattered photons is also presented. Finally an independent calculation of a simple 
warning flag for the occurrence of multiple scattering in a profile is described. 

Given the GLAS 532 nm channel specifications, Monte Carlo calculations show that the multiple 
scattering effect is expected to be significant (on the order of r|=0.6) in cloud situations, but is less 
than 20 percent (with rj=.83) in the worst aerosol situations. A method has been developed for 
calculations of GLAS multiple scattering using cloud and aerosol models based on a precise Monte 
Carlo radiative transfer model (Duda et al., 1999). However, a precise radiative transfer calculation 
to account for the effects of the multiple-scattered contribution is not practical for real time analysis, 
and approximation by semi-empirical methods is necessary. In any case, values of r| as a function 
of propagation depth and parameterized cloud and aerosol particle models can be calculated 
based on Monte Carlo calculations binned into look up tables for real time Level 2 processing. The 
defining equation used in this algorithm is 
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where F is the ratio of single scattering photons to multiple scattering photons referenced to a 
defined height and produced from Monte Carlo calculations. The defined height is usually the base 
of the particulate layer or the Earth’s surface, depending on the application. 


3.7.2. 1 Operational Multiple Scattering Correction Procedure 

The multiple scattering factor (rj) depends on the extent to which photons in the laser pulse have 
their trajectories altered by scattering events. This in turn is a function of the microphysical and 
physical properties of the cloud and aerosol layers in which the scattering occurs. Specifically, the 
effect of scattering depends on a) particle sizes within the scattering layer, b) the layer optical 
depth, c) the proximity of the scattering layer to the surface (for range-to-surface calculations only), 
and d) the physical thickness of the layer. It is important to understand that each of these factors is 
examined here independently, and the actual multiple scattering factor and scattering-induced 
range-to-surface delay are a result of both competing and additive affects from these various 
sources. 

Since the particle size for a layer cannot be directly observed or calculated for the GLAS data set, it 
is obtained through inputs into an effective radius look up table. The look up table is populated 
with seven specific particle sizes (0.6, 10.0, 12.0, 15.0, 22.0, 25.0, 40.0 microns) only. Size 0.6ja is 
reserved for aerosol layers and this is one of only two distinctions between clouds and aerosols in 
the multiple scattering correction process. The other particle sizes are based on clouds and are 
distributed globally, horizontal and vertical, on a monthly basis according to a pattern based on 
research by the GLAS lidar group (J. Spinhirne, personal communication). The multiple scattering 
table is computed from Monte Carlo runs at these effective radii only. The particle size table is 
dimensioned (18,18,12,3) where the indices represent latitude, longitude, month, and altitude zone 
respectively. 

Based on many Monte Carlo simulations at various optical depths, GLAS researchers have shown 
(Figure 3.7.3) that even though there is some small change of r| with increasing optical depth for a 
given layer physical thickness, using a constant value of r| for all OD is sufficient and will prevent 
erratic results. This is especially true since particle size is roughly estimated. For the Monte Carlo 
runs made to populate the multiple scattering look up tables, the layer optical depth was preset to 
0.47. 

Knowing the proximity of the layer to the Earth’s surface is needed for calculations of the multiple 
scattering factor when tied to the 1064 nm surface-reflectance-based column OD. The layer base 
height is binned in eight categories with no interpolation (0.10, 0.35, 0.75, 1.50, 3.00, 6.00, 10.00, 
14.00 km). Because the operational 532 nm lidar inversion processing is reference to the bottom 
of the layer, the height of the layer base is not used as an input into the 532 nm multiple scattering 
factor look up table. 

The physical thickness of the layer is an important input to Monte Carlo runs making up the 
multiple scattering factor look up table. Layer thickness is binned in seven categories with no 
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Figure 3.7.3 Plots of eta for 1064 nm as a function of optical depth based on Monte Carlo 
simulations. The particle size is set to 12 microns and the four graphs show layer thicknesses of 
0.1, 0.3, 0.8, and 1.5 km, respectively. The multiple curves in each graph mark various cloud base 
heights. 

interpolation (0.1, 0.35, 0.75, 1.50, 2.75, 4.75, 7.00 km). Interpolation would be problematic and 
probably not useful. Monte Carlo runs were made for each of these categories. 

Monte Carlo runs used to populate the look up tables executed a precise radiative transfer model 
based on a particle scattering phase functions simulating ice spheres (for clouds) and non-water 
particles (for aerosols) from Mie theory. T he ^ tables are dimensioned (7,7) for 532 nm with 
indices for particle size and layer thickness and (7,7,8) for 1064 nm with indices for particle size, 
layer thickness, and layer base height. Figure 3.7.4 displays various r\ results from the 532 nm 
look up table based on changing inputs. 
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Figure 3.7.4 Relationship of eta with particle size (x axis) and layer thickness (selected curves) for 
the operational 532 nm multiple scattering factor lookup table. 

3. 7. 2. 2 Operational Range Delay Calculation Procedure 
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The range-to-surface delay (range delay) is not used by the atmospheric processing routines, but 
is an important output on the GLA1 1 standard product for use by the altimetry group to properly 
correct the surface return height in the surface waveform algorithm from the effects of multiple 
scattering. The range delay is gotten from a table look up of calculations from many Monte Carlo 
simulations, much the same way as the multiple scattering factor table described in Section 
3.7.2. 1. The differences are that 1) the altimetry channel’s wavelength is 1064 nm, so the process 
is run for this wavelength, 2) the range delay is ground referenced rather than cloud base 
referenced, and 3) range delay is a linear function of optical depth, so the look up table is 
populated by optical depth conversion factors, which are slopes to the linear relationship. These 
slopes are calculated based on the differences between results from 0.45 and 1 .10 OD. 

The inputs to the range delay slope look up table for each layer are effective particle radius, layer 
thickness, and layer base height, identical to the multiple scattering table for 1064 nm The 
dimensions for the range delay slope table are (7,7,8) with indices for particle size, layer thickness, 
and layer base height respectively. To convert to range delay (5) in cm, use the following 
equation: 

(3.7.14) 5 = ar , 

where a is the range delay to optical depth slope conversion from the look up table and x is the 
layer optical depth calculated for the layer in the lidar inversion and corrected for multiple 
scattering. The full range delay for the atmospheric column is estimated by adding the range 
delays for each layer (both cloud and aerosol). Studies have shown that this technique compares 
favorably to calculating the full range delay directly. In order to be useful to the altimetry group, the 
range delay for the GLA1 1 product is converted to altimetry offset (mm) by multiplying by -10. The 
altimetry offset is set to invalid if the last (lowest) layer range delay is invalid. The offset is set to 
9999.9 if there is no lidar ground stroke detected. 

3.72.2 Multiple Scattering Warning Flag Calculation 

The multiple scattering warning flag is based on the 532 nm total column optical depth (aerosol 
plus cloud) calculated in the lidar inversion. It is intended as a way to quickly obtain information 
about the potential severity of multiple scattering with regards to the range-to-surface distance 
retrieval calculated by the altimetry processing software. It is output on the GLA1 1 product for use 
by the altimetry group. The range delay discussed in Section 3.1. 2.2 is a more rigorous and 
quantitative determination of the effect of multiple scattering on the range-to-surface distance. The 
multiple scattering warning flag has values ranging from 0-14, based on binning the total column 
optical depth into fifteen categories. The categories are presented in table 3.7.1 . The first 
category (0) is for total optical depth of 0.01 or less and the highest category (14) is for total optical 
depth of greater than 2.00. Category 1 5 is for invalid. 
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Table 3.7.1. Calculation of the multiple scattering 


Warning flag 


Optical Depth 

Multiple Scattering 
Warning Flag 

<0.01 

0 

0.01-0.03 

1 

0.03-0.06 

2 

0.06-0.10 

3 

0.10-0.15 

4 

0.15-0.225 

5 

0.225-0.30 

6 

0.30-0.40 

7 

0.40-0.50 

8 

0.50-0.67 

9 

0.67-0.90 

10 

0.90-1.20 

11 

1.20-1.60 

12 

1.60-2.00 

13 

>2.00 

14 


In summary, the algorithm will produce following quantities which will be written out to the GLA1 1 
data product: 

a. Multiple scattering factor (ranges from 0 to 1 ) 

b. Surface range delay estimate (millimeters) 

c. Surface range delay uncertainty estimate (millimeters) 

d. Multiple Scattering Effect Warning Flag (ranges from 0 to 1 5) 

e. Particle sizes estimated and used in the scattering calculation 

2.1 2 A Maximum Range-to-Surface Delay 

As mentioned above, the greatest uncertainty in determining the effects of multiple scattering 
arises from uncertain values of global particle size distributions. Using data from regional 
experiments, cloud and aerosol particle sizes are obtained in a fairly coarse grid, both spatially and 
temporally. In many regions of the world, such approximation is nonetheless the best that is 
possible. Particle sizes measured at a coastal Antarctic station, such as Palmer, for example, must 
be attributed to much of the coast of that continent. Broad latitudinal and longitudinal grids are 
defined within which the particle sizes are estimated from such regional studies. 

An alternate approach to understanding the effect of multiple scattering on the range-to-surface, 
however, exists. By this method, described below, it is possible to estimate the largest likely error 
in observations that will result from multiple scattering, even when particle sizes are unknown. In 
effect, this yields an upper limit to the uncertainty (due to multiple scattering) in the GLAS altimetry 
measurements. Though the approach is described here, the actual implementation will be done in 
the altimetry processing code. We will provide the altimetry team with a pre-calculated lookup table 
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which they will use to asses the maximum range-to-surface error. The following discussion is 
based on work performed by Ed Eloranta at the University of Wisconsin. 

Photon delays caused by scattering are dependent on both particle size and cloud altitude. 

Because particle sizes are not known, prudent estimation of ranging errors requires the assumption 
of worst case particle sizes at each altitude. The largest multiple scattering delays occur for particle 
sizes between 1 and 20 microns. Since particle sizes in this range are common in the atmosphere, 
the safest way to compute error bars describing the effects of multiple scattering would be to 
assume that clouds always contained particle sizes producing the maximum delay. In this way we 
are assured that this is a maximum upper limit on the ranging error due to multiple scattering. In 
addition, if the particle size at each altitude is selected to produce the largest delay, the delay is not 
a strong function of altitude, especially for small optical depths. 

Let us define the 1/e half width of the received surface return pulse (in the absence of atmospheric 
scattering) as: 8 p = Jsf+sf where So is the 1/e half width of the emitted laser pulse and 8 S is 
the standard deviation of the surface roughness elements. An increase in surface roughness 
produces effects which are indistinguishable from those caused by lengthening the laser pulse (So). 
As we know, atmospheric multiple scattering tends to broaden the received surface return 
waveform. Surface roughness and slope will have a similar effect. Surface roughness broadens 
both the directly reflected pulse and the multiply scattered return symmetrically and thus has no 
effect on the centroid estimate of altitude. However, because increased surface roughness 
produces a longer directly reflected pulse with a flatter peak, adding the temporal skewed 
atmospheric scattering will shift the center position of the Gaussian fit. 


Path delays computed from Gaussain fit, tau= 0.1 



Figure 3.7.5. Path delay (cm) computed from Gaussian fit with S p = 2 and t = 0.1 . 
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An examination of Figure 3.7.5 shows that while the particle radius that produces the worst delay 
changes with altitude, the maximum value of the photon delay is nearly constant with altitude. This 
simplifies the error calculation because the worst case error becomes a function of the total optical 
depth and it is not necessary to know the vertical profile of the scattering cross section. Using this 
observation, it is possible to generate a plot of the worst case errors as a function of surface 
roughness and optical depth which may then serve as a basis for assignment of realistic multiple 
scattering error bounds for the GLAS altitude measurements. 


Altitude error as function of optical depth and surface roughness 



Figure 3.7.6. Contour plot of altitude error as a function of optical depth 
and surface roughness. Note that the y axis is actually in units of meters, not cm. 

Figure 3.7.5 shows the maximum ranging error as a function of atmospheric optical depth and 
surface roughness. The figure provides a means of assigning an upper bound to the scattering 
induced error in the GLAS altimetry measurements. The optical depth is estimated from the 532 
nm lidar channel. The width of the received pulse determined from a Gaussian fit to the received 
waveform would be used as an upper limit on the value 8 P . The lidar science group will do the 
necessary calculation to prepare a two dimensional table that will have optical depth as one index 
and S p as the other as shown in table 3.7.2 below. The contents of the two dimensional array will 
be the maximum delay for each 8 P , optical depth pair. The values may be linearly interpolated from 
the given points to obtain the delay for specific values of optical depth and S p . We will supply this to 
the altimetry software group and they can use their measured surface pulse width and our total 
optical depth (output on GLA1 1) to index into the table to retrieve the maximum delay. 



Table 3.7.2 The calculation of the l,J indices of the maximum 


range-to-surface delay table 


6 P (m) 

1 Index 

Optical Depth 

J Index 

0.20 

1 

0.05 

1 

0.40 

2 

0.10 

2 

0.60 

3 

0.20 

3 

0.80 

4 

0.30 

4 

1.0 

5 

0.40 

5 

1.5 

6 

0.50 

6 

2.0 

7 

0.60 

7 

2.5 

8 

0.70 

8 

3.5 

9 

0.80 

9 

4.0 

10 

0.90 

10 

5.0 

11 

1.0 

11 

6.0 

12 

1.2 

12 

7.0 

13 

1.4 

13 

8.0 

14 

1.6 

14 

9.0 

15 

1.8 

15 

10.0 

16 

2.0 

16 
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4 Practical Application 


This section will address the practical issues related to coding, implementing and running 
the algorithms. These topics include the type and source of input data required to run 
the algorithm, execution time, program flow considerations (execution order), and 
examples of output where appropriate. Each algorithm will be addressed separately, 
and in the same order they were presented in section 3. 

4.1 Normalized Lidar Signal 

4. 1. 1 Required Input Data 

In addition to the raw lidar return signal for each channel, the normalized lidar signal 
(GLA02) algorithm will require the laser energy (reported at 40 Hz) and the two 
background measurements. Also required are the dead time correction table for the 532 
photon counting channel, the 1064 amplifier gain and attenuation settings, the 1064 
voltage offset and the factor relating digital counts to volts. While not explicitly used in 
the algorithm, the 532 channel etalon filter settings should be supplied, as it may be 
needed in subsequent processing (GLA07). The first data bin of the 532 channel is 
supposed to be 41 km above local DEM. In order to compute the range (R) used in 
equations 3.1 .1 and 3.1 .2, we need to know the range from the spacecraft to the top 
bin of the lidar profile. While not required, it is assumed that Global Positioning System 
(GPS) time and position (latitude and longitude) will be provided in the input data stream. 
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4.1.2 Algorithm Implementation 


This algorithm is relatively easy to implement and does not require a large amount of CPU time. 
There will be 3 main loops for the processing of the 3 distinct data layers (2 layers for the 1 064 
channel). The data in the upper two layers for the 532 channel must be normalized by the number 
of shots summed in that layer before being used in equations 3.1 .1 and 3.1 .2 (as will the top layer 
of the 1064 channel). The bottom layer of both channels do not need normalization as they are 
single (un-summed) shots. Note that the laser energy must be computed at 5 Hz for the 1 064 and 
532 channels in order to compute P’ from the 8 shot summed data (the 1 0 to 20 km layer). For the 
532 channel, the laser energy at 1 Hz must also be computed in order to use equation 3.1.1 in the 
upper layer. Similarly, the background must be computed at 5 and 1 Hz in order to process the 
data in the 10 to 20 km layer and the 20 to 40 km layer (532 only), respectively. These calculations 
should take place as the 40 Hz data are being processed, using the 40 Hz laser energy and 
background values. 

The saturation flag applies only to the 532 channel and will take the form of profiles, each bin of 
which will have a one-to-one correspondence with the data bins of the 532 channel. It will be a one 
byte value, where zero indicates that the 532 channel is not saturated and 1 denotes detector 
saturation. Whether a 532 bin is saturated will be determined from the magnitude of the 532 signal 
in that bin. The level at which saturation occurs has been determined during calibration procedures 
in the laboratory. 

The predicted height of first cloud top will be calculated at 5 Hz using the raw 532 channel data, 
which means that an 8 shot average of the data from the lowest layer must be computed. This 
would be tacked on to the corresponding normalized 8 shot sum form the middle layer to form a 
profile from -1 to 20 km (at 5 Hz). A search would then begin as described in section 3.1 to find the 
first bin which exceeds a cloud threshold level, which will be about 2-3 photons per bin (after 
background subtraction). The bin number that exceeded the threshold will be stored and the height 
will be calculated as described in section 3.1 .1 . The 5 Hz profile that is used to find the height of 
the first cloud top is not saved as output. It is discarded. 

As noted in section 3.1.1 .3, GLA02 will also be computing averages of the normalized signal at the 
calibration heights. A crucial part of this process (for the lower calibration height) will be cloud 
screening the data so that no shots that have clouds in them will be included in the average. The 
cloud search described above will be used for this purpose. If a cloud is found anywhere above 7 
km, then that second of data is not used for the computation of the signal average. 

After the normalized signal for the 3 layers is computed using 3.1 and 3.2, it must be scaled to fit in 
a 4 byte signed integer. A signed integer is required because P’ can be negative (due to the 
subtraction of the background value). The fairly large dynamic range of the computed signals 
warrants the use of a four byte integer. The scaling can be accomplished by applying a simple 
multiplicative scaling factor. 

4.1.3 Interpreting the Output 
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The output from GLA02 consists of calculated parameters as well as passed-through quantities 
which are not calculated or used by the algorithm, but which will be used in the creation of level 2 
data products. A list of all GLA02 output follows: 

1 . P at 532 nm for 3 layers: -1 to 1 0 km (40 Hz), 1 0 to 20 km (5 Hz), and 20 to 40 km (1 Hz) 

2. P at 1064 nmfor2 layers: -1 to 10 km (40 Hz) and 10 to 20 km (5 Hz) 

3. .532 nm saturation flag for the 3 layers: -1 to 10 km (40 Hz), 10 to 20 km (5 Hz), and 20 to 40 km 
(1 Hz) 

4. Range from spacecraft to top of data 

5. Predicted height of first cloud top (5 Hz) 

6. Ground return flag (bin number) and maximum ground signal 

7. 532 background (40, 5 and 1 Hz) and 1064 background (40 and 5 Hz) 

8. 532 laser energy (40, 5 and 1 Hz) 

9. 1064 laser energy (40 and 5 Hz) 

1 0. 532 laser energy quality flag (40 Hz) 

11.1064 laser energy quality flag (40 Hz) 

12. 532 integrated return from 40 to 20 km (1 Hz) 

1 3. 532 quality flag (1 Hz) - based on 1 1 above 

14. 1064 nm programmable attenuation setting (1 Hz) 

15. 532 nm etalon filter parameters (1 Hz) 

16. GPS time (1 Hz) 

Items 1 through 8 have been thoroughly discussed in section 3.1. The 532 and 1064 channel laser 
energy quality flags (9 and 10) are discussed in section 4.1 .5 below, as is the 532 integrated return 
and associated quality flag (11 and 12). The 1064 programmable attenuation setting (13) is used to 
compute the normalized signal (Equation 3.1 .2) and should be reported in the GLA02 input data 
stream. The 532 nm etalon filter parameters are as yet unspecified but should provide a measure 
of how well tuned the etalon was to the laser frequency. The etalon filter parameters are not used 
in GLA02, but may be needed in subsequent processing. 

4.1.4 Quality Control 

At this early stage of data processing quality control should be directed at assessing the health of 
the instrument and the fundamental soundness of the lidar return signal. The health of the laser 
can be assessed by monitoring the laser energy. For each channel, a quality flag is set for every 
shot (40 Hz) which characterizes the laser as follows: 

1 = full laser energy (within 90 percent of expected max value) 

2= marginal laser energy (between 90 and 70 percent of expected max value) 

3 = deficient laser energy (less than 70 percent of expected max energy) 

Used in conjunction with the above flags, the boresite can be assessed by integrating the 532 
return signal from 41 to 20 km altitude (l s ). This could be done using raw photon counts after the 
background (computed from equation 3.1 .3) is subtracted out. A quality flag should be set much 
like for the laser energy flag above, depending on the magnitude of the integrated return. Based on 
simulations, the expected number of integrated photons (the summation of 268 bins) per second 
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from this region of the atmosphere is about 1 900. For the 532 channel a quality flag can be 
formulated as follows: 

1= excellent signal strength (l s > 1800) 

2= good signal strength (1400 < l s < 1800) 

3= marginal signal strength (1000 < l s < 1400) 

4= poor signal strength (500 < l s < 1000) 

5= bad data (l s < 500) 

These limits are based on simulations and may be adjusted up or down after actual data are 
acquired. 

4.2 Attenuated Backscatter Cross Section 

4.2.1 Required Input Data 

This algorithm requires the output from GLA02, the normalized lidar signal and associated output 
as described in section 4.1 .3. In addition to this, the GLA07 algorithm requires MET data in order to 
compute the molecular backscatter cross section at the various calibration heights (see equations 

3.2.1 to 3.2.5). Realizing that MET data will most likely not be available at all times, it is important 
to also provide the standard model atmosphere as a backup source to obtain the required 
temperature and pressure profiles. The standard atmosphere actually consists of 5 models defined 
roughly by latitude and season as follows: Mid latitude and arctic for both summer and winter and 
tropical (5 standard atmosphere models in all). Additional required input includes the precision orbit 
determination (POD) data which includes latitude, longitude, spacecraft altitude and pointing angle. 
In order to transform the data height coordinate from above local DEM to above mean sea level, 
we need to know the spacecraft altitude (above mean sea level) and the range from the spacecraft 
to the top of the lidar profile. Also, to correct for vertical errors introduced by off-nadir pointing, the 
pointing angle is required. These issues are discussed fully in section 3.2. 1.1. 

4.2.2 Algorithm Implementation 

The main function of this algorithm is to compute and apply the lidar calibration constant to the data 
to form a continuous 5 Hz profile of attenuated backscatter cross section from 41 to -1 km for the 
532 channel, and from 20 to -1 km for the 1064 channel (the altitude is with respect to mean sea 
level). In addition, 40 Hz profiles from 10 to -1 km will also be generated for both channels. The 
calibration constant (C) will be computed from the signal average file generated by GLA02 and 
described in section 3.1 .1 .3. This is done immediately by GLA07, before anything else is done. 

This process will generate a calibration constant to be used for each second of the granule via 
interpolation between points or simply an average of all the calibration values. 

The ability to compute accurate C values for the 1 064 channel is in doubt. The algorithm will 
perform the calculations for the 1 064 channel as described in section 3.2.1 .2, but, at least initially, 
the laboratory calculated 1064 calibration constant will be used in equation 3.2.6. If subsequent 
analysis indicates that the 1064 C calculated from the flight data is good, then it may be used. 

Thus, for the 1064 channel, a flag should be built in which tells the software to use the laboratory C 
value or the C calculated from the atmospheric data. 
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The calculation of the lidar calibration constant requires the construction of accurate molecular 
backscatter profiles through the calibration layer(s). Since the entire (0 to 41km altitude) molecular 
backscatter profiles will be required by other GLAS atmospheric data product modules, it makes 
sense to compute them here. When this is done from MET data or from standard atmosphere data, 
it will be required to interpolate between the standard pressure levels to a vertical resolution 
equivalent to the lidar profile (76.8 m) as discussed in section 3.2.1 .2. Since the molecular 
scattering depends only on atmospheric density, it makes sense to first compute the density from 
the temperature, relative humidity and pressure at the geometric heights corresponding to the 
standard pressure levels and then use the hypsometric formula to compute the density between 
the standard heights. This will result in a smooth density profile with 76.8 meter vertical resolution. 

The output consists of 5 Hz full profiles (-1 to 41 km for 532 and -1 to 20 km for 1064) and 40 Hz 
profiles of only the lowest layer (-1 to 10km) for both channels. In order to form the 532 full profile, 

8 shots of the lowest layer are averaged and the corresponding 8 shot sum (after being 
normalized) of the middle layer (10 to 20 km) is placed above that, with the normalized 40 shot 
sum profile (20 to 41 km) above that. The same 20 to 41 km profile is used repeatedly for each of 
the 8 shots for a given second. The result is a 5 Hz full profile from -1 to 41 km with respect to 
mean sea level. For the 1 064 channel, 8 shots of the lowest layer are averaged and combined with 
the normalized 8 shot sum of the middle layer to form one 5 Hz profile from -1 to 20 km. Note also 
that the 532 saturation flag full profile (5 Hz, -1 to 41 km) must also be formed. In this case, instead 
of averaging in the lowest layer, we sum up the saturation flag for the eight shots yielding a number 
between 0 and 8. This is then combined with the saturation flag from the middle layer and the 
saturation flag profile from the upper layer (which is repeated 8 times). 

4.2.3 Interpreting the Output 

The output of GLA07 consists of profiles of calibrated attenuated backscatter cross section and the 
calibration constants for both channels. The 532 channel consists of 5 Hz, -1 to 41 km (548 bins) 
and 40 Hz, -1 to 10.3 km (148 bins) profiles which have had saturated bins replaced with estimated 
cross section provided by the 1064 channel (if a flag indicates that this should occur). The 1064 
channel output will consist of 5 Hz profiles from -1 to 20.5 km, and 40 Hz profiles from -1 to 10.3 
km. An example of one 5 Hz profile (532 nm) output from the algorithm is shown in Figure 4.2.1 . 
Examination of Figure 4.2.1 shows that the measured backscatter profile closely follows the 
molecular return until a thin cirrus cloud is encountered at 1 1 km. Below this layer, the signal 
continues to follow the molecular backscatter until the top of the PBL is reached (at about 1 .5 km 
height). Scattering within the PBL is much larger than the molecular scattering level due to the high 
concentration of aerosol there. Also seen in the figure is the large ground return signal which is at 
the 1 km height for this segment of data. 
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GLAS 532 nm Average Profile 



Figure 4.2.1 . An example of a 5 Hz 532 nm attenuated backscatter profile from 40 to -1 km which 
is output from the GLA07 algorithm. The thick black line is the attenuated molecular return at 532 
nm. 

In Figure 4.2.2, we have assembled many such profiles together and presented them in image 
form. This is probably the most informative way to display the data because it contains so much 
information. At a glance, one can see the various cloud layers, the boundary layer height and 
structure and any elevated aerosol layers that might be present. 



[3] 532 nm CAB (1/(m-sr)) 06 — Oct — 2003 2 1:13:04 — 22:49:35 GMT 



Calibration Constant: I.OOe+13 Latitude/Longitude Calculated Constant: 6.77e+12 


[4] 532 nm CAB (1/(m-sr)) 06 — Oct — 2003 22:49:36 — 23:59:59 GMT 



137 134 1 28 tIO —3-* --*6 —50 —5* -58 —63 O 

Latitude/ Lonqitude 

Calibration Constant: I.OOe+13 Calculated Constant: 6.37e+12 


Figure 4.2.2. 532 nm attenuated backscatter cross section displayed in color image form. The 
image is comprised of 2000 separate 5 Hz profiles like the one shown in Figure 4.2.1. 
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In addition to the backscatter profiles, the calculated calibration constants for both channels at the 
two heights are output as are the actual calibration values used in equations 3.2.12 and 3.2.13. 
Remember that the program is designed to use either the calculated value, a previously calculated 
C value or the laboratory calibration constant. With each calibration point a flag will be generated 
which characterizes the background condition (day, night or undetermined) which existed during 
the calculation of that C value. The 532 saturation flag profiles will be output at 5 Hz for — 1 to 41 
km and 40Hz for -1 to 1 0 km as described in section 4.2.2. 

The profiles of attenuated backscatter cross section, which are the main output from GLA07, will 
consist of 548 bins, each 76.8 meters wide and stretch from 41 km to -1 km above mean sea level. 
The process of shifting the bins to compensate for the varying topography will mean that some of 
the data will be cut off on top and some will be buffered with a missing data value at the end of the 
548 bin profile. For example, if data are being acquired over a region which is 5 km above mean 
sea level (as determined from the onboard DEM), the resulting acquired 532 nm profile will actually 
cover the region 46 to 4 km above mean sea level. The profile which will be output from GLA07 will 
truncate the 5 km of data above 41 km, and fill the region of the profile below 4 km with a missing 
data value (-999 is suggested). The same is true of the 1064 channel, except it extends to only 20 
km above mean sea level. Note that there will be a small percentage of time where the data are cut 
off at the bottom of the profile and padded at the top. This would occur when the DEM value is less 
than mean sea level. A complete list of the output for GLA07 follows: 

1 . 532 nm attenuated backscatter cross section, 41 to -1 km above mean sea level at 5 Hz 

2. 532 nm attenuated backscatter cross section, 10 to -1 km above mean sea level at 40 Hz 

3. 1064 nm attenuated backscatter cross section, 20 to -1 km above mean sea level at 5 Hz 

4. 1064 nm attenuated backscatter cross section, 10 to -1 km above mean sea level at 40 Hz 

5. 532 nm saturation flag profiles, 41 to -1 km above mean sea level at 5 Hz and 10 to -1 km at 
40 Hz. 

6. 532 nm calibration constants - upper (C30), lower (Ci) , and the actual C value that was used 
in equation 3.2.12 

7. 1 064 nm calibration constants - Ci , and the actual C value that was used in equation 3.2.1 3 

8. Calibration constant day/night flag (see discussion, section 3.2.1 .2) 

9. Calibration constant quality flag (see discussion section 4.2.5, below) 

10. Ground return bin as determined from POD and DEM 

1 1 . Predicted height of first cloud top, 5 Hz 

12. Ground return flag (bin number) and maximum ground signal 

13. 532 nm background at 40 and 5 Hz 

14. 1064 nm background at 40 and 5 Hz 

15. 532 laser energy at 40 and 5 Hz 

1 6. 1 064 laser energy at 40 and 5 Hz 

17. 532 laser energy quality flag at 40 Hz 

18. 1064 laser energy quality flag 

19. 532 nm integrated return from 40 to 20 km at 1 Hz 

20. 532 quality flag (1 Hz) based on 17 and 19 above 

21. 1064 programmable gain amplifier setting (1 Hz) 

22. 532 nm etalon filter parameters (1 Hz) 

23. GPS time (1 Hz) 

24. Precision Orbit Determination (POD) data (1 Hz) 
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25. Attitude information (pointing angle) 

26. Onboard Digital Elevation Model (DEM) value used (1 Hz) 

Items 1 through 10 are calculated by GLA07. The remaining output is either from the output of 
GLA02 or another input source. 

4.2.4 Quality Control 

Quality control should be implemented during the calculation of the calibration constant by 
checking the data quality flags that were generated by GLA02. Specifically, the laser energy flag 
and the integrated (20 to 40 km) return flag should be used to eliminate bad shots. Based on these 
flags, a bad shot counter should be kept during the calculation of C. If the number of bad shots 
exceeds say 5 percent of the total number of shots expected to be processed for a given 
calibration cycle (there are about 27,000 shots per 1/8 orbit), then the C value calculated for this 
calibration cycle should be flagged as questionable. Also time continuity must be checked during 
the calculation of C to check for large time gaps in the data that might adversely affect the 
calculation of the calibration constant. If a time gap greater than 30 or 40 seconds is encountered 
(total time for 1/8 orbit is about 675 seconds), the calibration constant should similarly be flagged 
as questionable. 

The quality of the calibration constant can be assessed by looking at its variability with time and the 
difference between the constants calculated at the two different heights. The 532 nm attenuated 
backscatter cross section profiles can be checked by normalizing them by the attenuated molecular 
profile. This should produce a profile that ranges between 0.9 and about 10.0. This test could only 
be applied to data with a ground return as the values below thick clouds will approach zero. 

A test that could be applied to all the data would be to integrate the attenuated backscatter ((3’) 
from 40 to 20 km and divide by the integrated attenuated molecular backscatter (p m T m 2 ) to form a 
ratio that should be very close to unity. A major deviation from one would indicate a problem. 

4.3 Particle Layer Height and Earth's Surface Height (GLA09) 

The implementation of the algorithms to find vertical cloud and aerosol layer boundaries and the 
height of the earth's surface (ground height) will require only modest resources in terms of coding 
and execution time. The processing will be done on a time series of 4-second segments. Each of 
these will be composed of matrices of values consisting of on the order of 1 0000 points. Only 
elementary arithmetic computational and logical functions and testing will be done on these to 
produce the output. The output consists of cloud layer height at 4 second, 1 second, 5 Hz and 40 
Hz for a maximum of 10 layers at all 4 horizontal resolutions. Also generated are elevated aerosol 
layer heights at 4 second resolution (maximum of 8 layers) which are stored on the GLA08 product. 

4.3.1 Required Input Data 

The vertical boundaries of the horizontal surfaces of particulate layers and the earth's surface will 
be found by testing profiles of attenuated backscatter coefficients (backscatter cross-section) 
against thresholds developed from the profiles themselves. The profiles will be those developed in 
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GLA07. Both 532 nm and 1064 attenuated backscatter profiles are used. The following list 
summarizes the required input for each 4-second granule. 

a) 160 40 Hz attenuated backscatter coefficient profiles, -1-10 km, from GLA07; 

b) 20 5Hz attenuated backscatter coefficient profiles, -1-40km, 

c) 4 sets of 1 -second DEM values corresponding to the 4-second processing interval, with each set 
containing the mean, maximum and minimum values of the height of the earth's surface in the local 
1 -degree square grid 

d) one current atmospheric profile, 0-22km, containing pressure, temperature, height, and relative 
humidity. 

4.3.2 Algorithm Implementation 

Cloud and aerosol layer boundary searches will be limited to the lowest 22 km above ground level. 
For each 4-second interval, the computations proceed in the following manner. 

The input profiles are acquired. These consist of 160 40Hz profiles and 20 5Hz profiles. Each of 
the 40Hz profiles extend from -1 to 10 km. The total number of samples therefore is 
{[1 0-(-1 )] km/ 0.0768 km}= 143 samples/profile; 

143 samples/profile x 160 profiles=22880 samples. 

Each 5 Hz profile extends from -1 to 40 km; thus, each has 533 samples. So, the total number of 5 
Hz samples is 

( 533 samples/profile) x 20 profiles = 10660 samples. The total number of backscatter coefficient 
samples required for each 4-second interval is 22880+10660=33540. This value represents the 
most significant demand on computer memory storage for this algorithm. 

The first task will be to detect and position the signal from the earth’s surface (ground signal). As 
mentioned in Section 3.3.1 .2, the ground height for the 1 Hz and 0.25 Hz profiles will be the 
average of the ground height from the appropriate 5 Hz profiles. The details of the ground 
procedure are found in that section. Since the algorithm requires that cloud boundaries be found at 
the 4-second resolution first, a 4-second average profile will be found by averaging the 5 Hz 
profiles. The averaging will extend from -1 km to 22 km. When a layer is found at the 4 second 
resolution, it will be then be characterized as either cloud or aerosol using a set of discrimination 
criteria described in section 3.3.1 .2. If it is aerosol, no further searches are done at the higher 
resolutions. If it is determined to be a cloud, the one-second profiles will be analyzed after the 
processing of the 4-second profiles is complete. 

The 4-second average profile will be divided into n s segments. These will not necessarily be of 
uniform length. In each of these segments, a minimum value will be found. A measure of the 
random noise associated with the molecular signal will be computed. For observations taken in 
sunlight, the standard deviation of the background signal will be used. This value will be found by 
using values in the final 1 3 samples of the profile, which occur after the laser pulse is extinguished 
by the ground. The lack of a significant background signal in night observations requires that the 
variability be estimated from 18-19 km. portion of the profile. The atmosphere at that altitude will be 
free from strong non-molecular scattering constituents. Therefore, the variability of the signal found 
there will be representative of the variability found in the particulate free portions of the entire 
profile. A constant fraction of this variability will be added to each of the segment minimums. The 
optimum value for the fraction will be determined from modeling studies. A threshold profile, 
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extending from 0-22 km, will be formed by linear interpolation and extrapolation of the segment 
points. 

Layer boundaries can now be found from the 4-second profile. Starting at the first sample, which 
represents the highest point, each profile value will be compared to the threshold profile. The 
presence of layer will be designated false. When a certain number of consecutive samples are 
found to be greater in value than the threshold, a top-of-layer will be located where the first sample 
exceeded the threshold. The presence of layer will be designated true. The comparison of profile 
values will continue downward. When certain number of consecutive values are found to be less 
than the threshold, a bottom-of-layer will be located where the first of the consecutive samples was 
found and the presence of layer will return to false. This process will continue to ground level. The 
location of the ground level will be the average of the ground levels derived from the 5Hz profiles, 
as described in Section 3.3.1 .8. If no ground signal was detected, the layer search algorithm will 
continue to the minimum of the 4 DEM minimums associated with the 4-second segment. When 
this procedure is complete, the set of layer boundaries and ground location of the 4-second profile 
will be known and stored. The result of the 4 second search will then be subject to the 
cloud/aerosol discrimination routine. If the layer is deemed an aerosol, further higher resolution 
processing is not done and there resulting layer height at 4 seconds is stored on the GLA08 
product. Up to 5 aerosol layers are stored below 20 km altitude and 3 above 20 km. 

If the 4 second layer is determined to be a cloud, then the procedure to find cloud boundaries for 
each of the 1 second profiles contained in a 4 second segment will be applied to the GLAS 
backscatter coefficient profiles. The algorithm will be the same as that described for the 4-second 
profile, with the following alteration. No boundaries outside of those defined by the results from the 
four-second profile (plus a small delta) will be accepted. If such a boundary or layer is found, it will 
be considered a false positive result caused by relatively larger random noise. Layers that are 
wholly outside of 4-second layers will be eliminated. Each 1 Hz result will be assigned a ground 
height computed from the average of the corresponding 5 Hz profiles. In a like manner, the layer 
boundary processing of the 5 Hz data will be same as the 0.25 Hz and 1 Hz. The 5 Hz layers will 
be required to exist in layers defined by the 1 Hz results. 

Finally, the 40 Hz profiles will be processed. Two factors force the processing to be somewhat 
different than that at the other frequencies. First, the 40 Hz data extend only from -1km to 10 km. 
Second, the low signal to noise precludes reliable detection of rarefied, optically thin clouds. In 
general, it is expected that only dense clouds will be reliably detected at 40 Hz. But knowledge of 
the location of cloud layer boundaries of these types of clouds in the lowest part of the atmosphere 
is very important for certain types of studies. For these reasons, the processing of 40 Hz data will 
proceed as follows. A ground height search algorithm will be applied independently to each of the 
profiles. The random noise factor will be estimated from the background portion of each profile. A 
threshold profile will be developed for the region 0-4km. If any cloud layers are detected, only the 
lowest one of those, confined to layers detected at 5 Hz, will be designated as a layer. 

4.3.3 Interpreting the Output 

The output of GLA09 will consist of the following products, for each 4 second processing segment: 
1 ) Results at 0.25 Hz frequency, 1 set 
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a) Vertical locations of the top and bottom of up to ten cloud layers, 0-22 km derived from the 
532 channel 

b) Vertical locations of the top and bottom of up to ten cloud layers, 0-20 km, derived from the 
1064 channel 

c) Quality flags for each layer 

d) Temperature, pressure and relative humidity associated with the top and bottom of each 
layer 

e) Ground height which will be the average of 20 5 Hz ground height results or indication of 
negative results if no ground was detected in the 4 second interval; 

f) Time, location and DEM information 

2) Results at 1 .0 Hz frequency, 4 sets 

a) Vertical locations of the top and bottom of up to ten cloud layers, 0-22 km. Derived from 
the 532 channel and confined to the layer boundaries detected at 0.25 Hz; 

b) Vertical locations of the top and bottom of up to ten cloud layers, 0-20 km, derived from the 
1064 channel and confined to the layer boundaries detected at 0.25 Hz 

c) Quality flags for each layer 

d) Temperature, pressure and relative humidity associated with the top and bottom of each 
layer 

e) Ground height which will be the average of 5 Hz ground height results or indication of 
negative results if no ground was detected in the 1 second interval; 

f) Ground height which will be the average of 5 Hz ground height results or indication of 
negative results if no ground was detected in the 1 second interval; 

g) Time, location and DEM information 

3) Results at 5 Hz frequency, 20 sets 

a) Vertical locations of the top and bottom of up to ten cloud layers, 0-22 km. Derived from 
the 532 channel and confined to the layer boundaries detected at 1 Hz; 

b) Quality flags for each layer 

c) Temperature, pressure and relative humidity associated with the top and bottom of each 
layer 

d) Ground height or negative results if no ground was detected; 

e) Time, location and DEM information 

4) Results at 40 Hz frequency, 1 60 sets 

a) Vertical locations of the top and bottom of one cloud layer, in the range 0-4 km, the lowest 
of any detected and confined to layer boundaries detected at 0.25 Hz (532 channel) 

b) Vertical locations of the top (no bottom) of one cloud layer, in the range 0-1 0 km, the 
highest detected (1064 channel) 

c) Ground height or negative results if no ground was detected; 

d) time and location information 

The tops and bottoms of layers are the heights h (in km above mean sea level) at which the layer 
signal becomes distinguishable from the molecular signal. In general, within the meaning of layer 
boundary at any of the time resolutions, the actual cloud boundary, h a , will be within a range of h- 
0.1 16km<h a <h+0.1 16km. 

If a ground signal is detected, than all layer boundaries are considered valid within the uncertainty 
limits. If no ground signal is detected, then the value of the bottom of the lowest layer has no 
meaning other than to indicate the height at which random noise first conceals the atmospheric 
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signal. Any layer of sufficient density and optical depth to cause multiple scattering to obliterate the 
location of the bottom of the layer will be assumed to fully attenuate the laser pulse. The bottom of 
the layer would not be meaningful in any such case. 

4.3.4 Quality Control 

The quality of the results of the GLA09 boundary procedure will be judged by how successful it is 
at finding all detectable cloud layers and locating their boundaries in the atmospheric profile. A 
significant advantage to the algorithm is that its application to a given time segment is independent 
of any GLAS observations outside of the segment. Quality of the results will be primarily controlled 
by the signal to noise ratio at any point in the profile. The results for each layer are superimposed 
on backscatter images and periodically reviewed. All results are assigned with quality flags based 
upon the ratio of the signal within the layer to that outside (above) the layer. 

The best way to judge the general quality of the results of the boundary algorithm is to plot the 
computed cloud boundaries on top of image segments constructed from lidar profiles. Such images 
reveal readily systematic and random faults in the results of the procedure. These inspections will 
be done on samplings of the results on a routine basis. If these reveal significant shortcomings in 
the method, the parameters used in the computation of thresholds will be adjusted to fix the 
discrepancies. 

4.4 PBL and Elevated Aerosol Layer Height (GLA08) 

4.4.1 Required Input Data 

The algorithm requires the 5 Hz profiles of the 532 nm attenuated backscatter and selected other 
components of the GLA07 output. These include the ground bin and various data quality flags. In 
addition, GLA08 requires the 5 Hz (high resolution) cloud boundaries output from the cloud 
detection algorithm (GLA09) and the 4 second aerosol layer heights (which are found by the 
GLA09 algorithm, but not output to the GLA09 product). Also required from GLA09 is the 5 Hz 
ground detection flag. Profiles of molecular backscatter cross section are also required since they 
are used to determine the bottom threshold as discussed in section 3.4.1 .2. 

4.4.2 Algorithm Implementation 

The algorithm can be implemented on any standard workstation with sufficient memory and CPU 
resources. To be most efficient, the 150 km record (100, 5 Hz profiles) of lidar data used to find 
the average PBL height should be kept in memory. The total memory required is less than 1 Mb, 
and the CPU requirements fairly minimal. Even though there are similarities between the PBL and 
particle layer algorithms, we have implemented them separately. For instance, the need for 20 
second and 4 second averaged profiles is common for both algorithms. However, the criteria for 
the composition of the 20 second averages is not quite the same. For the elevated aerosol layer, 
all shots are used regardless of the presence of a ground return or clouds. The PBL height, on the 
other hand, eliminates all shots without a ground return that have clouds above 6 km. 

4.4.3 Interpreting the Output 
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The output stored on GLA08 consists of planetary boundary layer height at high resolution (5 Hz or 
1 .5 km) and low resolution (0.25 Hz or 28 km). It will also contain the top and bottom height of a 
maximum of five elevated aerosol layers below 20 km at 0.25 Hz resolution, and a maximum of 
three aerosol layers above 20 km at a horizontal resolution of 0.05 Hz (150 km). When an aerosol 
layer is found above 1 0 km, and the temperature at the height of the layer is below -80 °C, and the 
latitude is poleward of 65 degrees, a Polar Stratospheric Cloud (PSC) flag is set to indicate a very 
high likelihood of the layer being a PSC. If the layer temperature is above -80 °C, but less than - 
70 °C, the flag is set to a different value to indicate a lesser likelihood of it being a PSC. The PSC 
flag will have the value of zero at all other times. All heights generated will be in kilometers above 
mean sea level. An elevated aerosol layer is defined as a region of increased lidar backscatter 
(above local ambient values) which has a minimum thickness of 230 meters (3 data bins). The data 
input to GLA08 will already have been processed by the cloud and aerosol height detection 
algorithm (GLA09) 

A partial list of the output of GLA08 follows: 

1 Aerosol layer top and bottom height above 20 km (max 3 layers) at 20 second resolution 

2 Aerosol layer top and bottom height below 20 km (max 5 layers) at 4 second resolution 

3 Aerosol layer height quality flags (4 and 20 seconds) 

4 Flag indicating whether given layer is a Polar Stratospheric Cloud (PSC) 

5 Planetary Boundary Layer (PBL) height at 5 Hz and 4 second resolution 

6 PBL height quality flags (4 seconds and 5 Hz) 

7 Clear/Cloudy flag for PBL height at 5 Hz and 4 second resolution 

8 Ground return bin (5 Hz and 4 seconds) 

9 Flag indicating whether GLA08 or GLA09 was used to produce the 4 second aerosol layer 
top and bottom heights 

10 Precision Orbit Determination (POD) data (1 Hz) 

11 GPS time (1Hz) 

12 Orbit Number 

13 PAD Pointing Vector (1 Hz) 

4.4.4 Quality Control 

Validation of the algorithm output can best be accomplished by overlaying the PBL and aerosol 
layer heights on top of the images, much like what is shown above. We have found from 
experience that visual inspection can reliably distinguish aerosol boundaries when the data are 
presented in this form. The lidar retrieved PBL heights can then easily be compared with the visual 
estimation of PBL height. The same is true for the elevated aerosol layers. Other validation 
approaches include using nearby radiosonde data to determine PBL depth and checking it against 
the lidar measurement (provided it is within a certain distance to the radiosonde station). Over land, 
it may be possible to use the MET data which is ingested by the GLAS ground processing system, 
and ground based lidar systems. 
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4.5 Optical Properties of Cloud and Aerosol Layers 
4.5.2 Required Input Data 

The algorithm that produces the GLA1 0 and GLA1 1 level 2 standard products will have as its 
starting point the 5 hertz GLA07 532 nm backscatter output. The operational processing of the 
atmospheric products is done in 20-second blocks. There are several parameters processed from 
the lidar signal or pulled from ancillary data sets that are used as inputs to produce extinction and 
corrected backscatter profiles and optical depth. 

4.5.1. 1 Retrieved Parameters from GLAS Lidar Signal 

1 ) 100 5-hertz average calibrated attenuated backscatter profiles stored in 20-second buffer 
containing the full vertical extent of the GLAS lidar data (41 km to -1 km). These profiles are 
then further processed to 1 -second, 4-second, and 20-second average profiles. Standard 
deviation profiles are calculated as follows: 1 -second averages use 5 hertz profiles, 4- 
second averages use 1 -second profiles, and 20-second averages use 4-second profiles. 

2) 5 4-second average cloud layer top height arrays (up to 1 0 layers vertically) from standard 
product GLA09 stored in 20-second buffer 

3) 5 4-second average cloud layer bottom height arrays (up to 10 layers vertically) from 
standard product GLA09 stored in 20-second buffer 

4) 1 20-second average stratospheric layer top height array (up to 3 layers above 20 km 
accompanied by a polar stratospheric cloud (PSC) flag) from standard product GLA08 

5) 1 20-second average stratospheric layer bottom height array (up to 3 layers above 20 km) 
from standard product GLA08 

6) 5 4-second average tropospheric elevated aerosol layer top height arrays (up to 5 layers at 
or below 20 km accompanied by a PSC flag) from standard product GLA08 stored in 20- 
second buffer 

7) 5 4-second average tropospheric elevated aerosol layer bottom height arrays (up to 5 
layers at or below 20 km) from standard product GLA08 stored in 20-second buffer 

8) 5 4-second average PBL top heights (up to 1 layer) from standard product GLA08 stored in 
20-second buffer 

9) 5 4-second lidar ground stroke heights from standard product GLA08 stored in 20-second 
buffer (height used is highest height in 4-second span) 

10) 100 5-hertz profiles of detector saturation status stored in a 20-second buffer. The 
saturation profiles are further processed into 1 -second, 4-second, and 20-second 
averages. 

4. 5. 1.2 Retrieved Parameters from Ancillary Data 

1) 20 1 -second average meteorological profiles of temperature, wind, relative humidity, and 
pressure from most current NCEP 1-km gridded model initialization matching ground track 
of satellite stored in 20-second buffer. The profiles are further processed to 4-second 
averages and 20-second averages. 
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2) 20 1 -second average 532 nm molecular backscatter profiles matching ground track of 
satellite (calculated from met profiles) stored in 20-second buffer. The profiles are further 
processed to 4-second averages and a 20-second average. 

3) Tropopause height based on latitude and season of year from standard atmosphere lookup 
table (1 value for 20 -second block) 

4) 20 1 -second average ground track latitude and longitude values stored in 20-second 
buffer. The averages are further processed to 4-second and 20-second averages. 

5) 20 1 -second average Earth surface height values based on 1-km Digital Elevation Map 
(DEM) stored in 20-second buffer. The averages are further processed to 4-second and 
20 -second averages. 

6 ) 20 1 -second averages of the laser tilt angle from nadir. The averages are further 
processed to 4-second averages and a 20-second average. 

7) Aerosol type assignments for each aerosol layer in 20-second buffer (The initial default 
separates out elevated layers from the PBL, PSCs from non-PSCs, and stratospheric from 
tropospheric. The sources of the type assignment for elevated tropospheric and the PBL 
are static latitude-longitude indexed global look-up tables of aerosol types based on 
climatological locations of aerosol source regions and transport preferences. For laser 
period L2A, a more sophisticated assignment is used based on the GEOS-4 global aerosol 
model initialization every 12 hours, separating out tropospheric and PBL assignments 
based on the vertical distribution of aerosols in the model and the GLAS PBL height.) 

The multi-scattering factor, 77 , (relationships formulated from section 3.7) will all be calculated 
based in whole or in part on pre-defined look up tables distinguishing between cloud and aerosol 
regimes. Elevated layers will be assessed for the capability of calculating S ' p from the signal 
profile. For those particulate layers where S' p can not be calculated, work done by Ackermann 
(1998) showed that reasonable estimates of S p for aerosols can be matrixed using location 
information (continental, maritime, and desert). Similar estimates can be done for clouds involving 
cloud phase and temperature. PSCs will be obtained from a subset of the aerosol matrix. S ' p can 
than be estimated by applying the estimate of 77 to the S p value. The following two sections 
describe the current default decision matrices of the S p look up tables in detail. 

4. 5.1. 3 Aerosol Extinction to Backscatter Ratio (S p ) Assignments 

Simple backscatter lidar inversion algorithms have to input an estimate of the extinction-to- 
backscatter (lidar or S) ratio in order to solve for optical properties. For every layer detected in the 
20-second block, a default S ratio is assigned. For some elevated layers where the S ratio can be 
calculated from the signal loss through the layer, the calculated S ratio is used. The algorithm to 
calculate the S ratio from signal loss is described in section 3.5.1 .1 . An S ratio use flag is used to 
keep track of the source of the final S ratio used. The ultimate accuracy of the optical properties 
solution will depend on the accuracy of the input S ratios. Because of the usually weaker signals of 
aerosol layers compared to clouds or the proximity of the layer to the ground, most aerosol layers 
are forced to use the default value in optical properties retrievals. Table 4.5.1 shows the cross- 
mapping of the GLAS S ratio assignments with aerosol type (see Item 7 under Section 4.5.1 .2) for 
elevated tropospheric and PBL layers. These values represent the latest in S ratio studies (as a 
function of aerosol type) in the literature plus case study results from AERONET, MPLNET, and 
CPL. It should be noted that these default S ratios are the “true” S ratios that would be used if 
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there were no multiple scattering (single scattering conditions). S ratios calculated from the signal 
loss through the layer are the “effective” S ratio, denoted by a hyphen (S’). A true S ratio must be 


Table 4.5.1 S Ratio Estimates based on Aerosol Type 


Aerosol 

Index 

Aerosol Type 

S-ratio 

Assignment 

0 

Default 

35.0 

1 

Sulfate+carbon 

67.5 

2 

Carbon 

62.0 

3 

Salt+dust 

32.5 

4 

Salt 

28.5 

5 

Sulfate 

60.0 

6 

Dust+carbon 

58.1 

7 

Salt+dust+sulfate 

47.2 

8 

Salt+carbon 

49.1 

9 

Salt+sulfate 

47.9 

10 

Dust 

42.5 

11 

Salt+dust+carbon 

48.2 

12 

Dust+sulfate 

56.5 

13 

Salt+carbon+sulfate 

53.3 

14 

Dust+carbon+sulfate 

58.9 

15 

All 

52.3 


converted to the effective S ratio when used inside the optical inversion algorithm. The simple 
relationship between real S and effective S is: 

(4.5.1) S\ = rjS p , 


where S p is the true S ratio of the particulate layer and □ is the multiple scattering fact 

Calculation of the multiple scattering factor of a layer is covered in a Section 3.7. The S ratio 
reported in the GLAS product GLA10 is the true S. 

PSC layers (both above and below 20 km) and non-PSC layers in the stratosphere have a different 
method for developing the default S ratio based on equations using backscatter strength and 
relative humidity. The S p ratio of both PSCs and regular stratospheric aerosol layers are based on 
backscatter strength (Gobbi, 1995) with different coefficients for Type I, Type II, and non-PSC 
stratospheric aerosol as shown in Eq. 4.5.2 

(4-5.2) , 


where p p is the aerosol backscatter cross-section in 1/cm-sr. This backscatter value can be 
estimated by finding the average value of (1 ,025*total attenuated backscatter - molecular 
backscatter) over the layer. In the GLAS algorithm, Type I is a PSC where the average relative 
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humidity is less than or equal to 95% and Type II is greater than 95%. Figure 4.5.1 shows the 
decision matrix used for stratospheric aerosols and the equation results of Type I and II PSCs and 


GLAS Stratospheric Aerosol Extinction to Backscatter Ratio (S p ) Default Matrix 



Figure 4.5.1 Flow diagram for GLAS Stratospheric S ratio Default Matrix 

non-PSC stratospheric aerosols. The tropopause height described in item 3 of Section 4.5.1 .2 is 
used to determine whether or not a layer is in the stratosphere by comparing to the layer top 
height. 

To separate Tropospheric and Stratospheric Layers: 

• Decide which standard atmosphere to use based on time of year and latitude (use 30N 
and 30S and Arctic Circle and Antarctic Circle for latitude boundaries plus use mid 
October-mid April and mid April-mid October as season boundaries). 

• T ropopause height (m)= 1 7000 for tropical 

9000 for arctic winter 
10000 for arctic summer 
14000 for mid latitude summer 
13000 for mid latitude winter 

4. 5. 1.4 Cloud Extinction to Backscatter Ratio (S p ) Assignments 

Cloud layers will be assigned a best default value of S p based on the matrix in Figure 4.5.2. The S p 
function shown in the figure is dependent on mean cloud layer temperature (degrees Centigrade) 
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and was created in-house by the GLAS science team based on available information from previous 
lidar results and ice crystal size and shape studies. For clouds whose temperatures are above -13 
C, the S ratio is 17.8 sr. 


GLAS Cloud Extinction to Backscatter Ratio (S p ) Default Matrix 



Average Cloud Temperature (C) 


Figure 4.5.2 Flow diagram of default cloud S p ratio assignments for use in optical property 
calculations. 


4.5.2 Algorithm Implementation 

The overall scheme of the aerosol optical properties procedure is to retrieve optical properties layer 
by layer sequentially from 41 km height to the surface at an along track resolution of 4-seconds (28 
km). This is then followed by the cloud procedure is to retrieve optical properties layer by layer 
sequentially from 20 km height to the surface at a resolution of 1 -second (7 km). In addition, in 
order to catch weak aerosol loading in the troposphere outside distinct layers, a procedure has 
been added to calculate the extinction profile and the cumulative optical depth of the “free 
troposphere” in the vertical zone between 20 km and the top of the first cloud or PBL, whichever 
occurs first. There are a maximum of nine potential aerosol layers in the column. The first three 
are reserved for layers found above 20 km. Layers 4-8 are reserved for elevated tropospheric 
layers. Layer 9 is reserved for any PBL detected. Layers are sequential inside their reserved 
space. In other words, if only an elevated tropospheric layer and the PBL were detected, then only 
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layers 4 and 9 would contain valid data. In order to raise signal-to-noise to an acceptable level, the 
optical retrievals above 20 km are calculated using 20-second averages, then this average is 
parceled out to the 5-four second retrievals in the 20-second block. In calculating the attenuation 
effect of the layers above on the current layer being analyzed, 4-second average cloud layers are 
optically analyzed in the mix, but only for their attenuation. There are a maximum of ten potential 
cloud layers. 

The end result of a 20-second output package should have 20 1 -second cloud backscatter and 
extinction profiles that contain calculated data in only those bins where cloud layers were found. 

All other bins contain defaults. The output package should contain 5 4-second aerosol backscatter 
and extinction profiles that contain calculated data in only those bins where aerosol layers were 
found. All other bins contain defaults. Layer quantities saved in output (clouds and aerosols 
separate) include Sp calculated from Integrated Ratio Technique (Sp=S’p/eta), Sp pulled directly 
from defaults, and a flag stipulating which Sp was actually used in the optical inversion, plus layer 
true optical depth 

A flow chart showing an overview of the optical parameter calculations is found in Figure 4.5.3. 

The critical component of the algorithm is the evaluation of the integral to compute y (see equation 
3.5.20). The flow of the algorithm proceeds as follows. For each profile P n (first the four second 
and then the one second resolutions), the levels where aerosol and cloud boundaries exist are 
obtained and differentiated. The molecular transmission to the top of the highest layer is computed 
to the 2Xsec0 power and used as k(zt) in equation 3.5.10. S p for the layer is computed (see 
section 3.5.1 .1) when the backscatter profile for a given layer is found to be appropriate for 
independent S ' p analysis. Otherwise a default value is used based on whether it is cloud, PSC, or 
aerosol. Obviously, whenever the S ratios are assigned rather than calculated, they could be 
higher or lower than the actual values. Miss-assigned values occasionally result in overshooting 
transmission thresholds during the processing of Equation 4.7. Too high of S ratio can be 
monitored because the calculated transmission is more likely to reach the minimum transmission 
threshold before processing reaches the bottom of the layer. Too low of S ratio is harder to get a 
handle on because it tends to decrease the rate of transmission reduction through the layer. The 
only scenario where it shows itself is the unstable case where it causes the calculated transmission 
to rise as you go deeper into the layer. The GLAS algorithm has added a feature where if the 
minimum transmission threshold is reached prematurely before processing has reached the layer 
base, then S' is reduced by 0.50 sr and processing of the layer starts over. This iterative process 

continues until processing reaches the layer bottom or the number of iterations exceeds 30. 
Knowledge of whether or not this iteration step occurred is kept in the S ratio use flag (1 =lookup 
default, 2=calculated, 3=modified default, 4=modified calculated). 

The integral for forward inversion is evaluated using a straight-forward rectangular summation. 

The terms of the summation are T^ x ~ l)secd P n Sz . The value of T p 2sec0 is computed for each 

level z in the layer. Computation for any subsequent layer will use the same method except that 
the k(zt) value will be re-computed as : 

(4.5.3) I B {z t ) = T' p 2 ™\z a )T 2 m x ™\z t ), 


109 



DFD 4.3: Calculate Level 2 Cross Sections and Optical Depths 



Figure 4.5.3 Flow diagram of level 2 optical parameter calculations. 


where T' 2sec0 (z a ) is the particulate slant transmission squared a the bottom of the layer above 
and T 2 (z t ) is the molecular transmission squared calculated down to the level of zt, the top 
location of the current layer. The backward inversion method initializes particulate backscatter to 
zero one bin below layer bottom and calculates corrected particulate backscatter, (3 p (z), at each bin 
based in part on the particulate and molecular backscatter at one bin below, starting at the layer 
bottom and ending at the layer top. Once the backscatter profile is calculated for the layer, 
extinction and transmission profiles can be calculated. The computer equations governing the 
backward inversion algorithm for each bin location (ib) in the layer are as follows: 

(4.5.4) A = (S' p - S m XP m ( ib ) + p m (ib + 1))AZ 


(4.5.5) p p (ib) 


P n (ib) * exp(Z) 


P„(ib + 1) 

(3 p (ib + \) + P m (ib + \) 


+ S' p (P n (ib + \) + P n (ib)^ V (A))^Z 


P m (ib) 


(4.5.6) a (ib) = ^j3 (ib) 

n 


no 


(4.5.7) 


T' p \ib) = 


p K m 

T~i(ib)(p p (ib) + p Jib)) 


The backward or forward inversion continues throughout each particulate layer as per the eight 
rules outlined above until T’J ce (z) < T r or the signal from the earth’s surface is detected. 


The algorithm is moderately computationally intensive. Results indicate that to process an orbit of 
data for the GLA1 0 and GLA1 1 products would take about 1 .05 minutes of cpu time. 

4.5.2. 1 Specific Optical Properties Retrieval Procedure for the Free Troposphere 

A feature that was added late to the GLAS optical processing algorithm involves the free 
troposphere. This is an attempt to calculate the background extinction in the free troposphere plus 
any unidentified aerosol missed beyond identified layers. The procedure is very similar to the 
generic algorithm discussed above except for the following: 

1 ) The top of region is set to 20 kilometers. 

2) The bottom is set to just above the top of the highest cloud, the PBL, or Earth’s surface, 
whichever is highest. 

3) The S ratio assignment always defaults to the elevated tropospheric aerosol index value 
for the current time and location. 

4) A constant custom value of the multiple scattering factor is set to 0.90. 

5) The processing is only done at night. 

No consideration for detector saturation is applied since the lidar signals retrieved are assumed to 
be well below saturation thresholds. The output of the free troposphere processing is limited to 
profiles of extinction and backscatter cross section, total aerosol optical depth of the free 
troposphere region as defined above, and the bottom height of the free troposphere region. Note 
the definition of the free troposphere region may include embedded elevated aerosol layers that 
have been separately identified and will be processed later in the layer by layer processing. The 
free troposphere aerosol optical depth will then be refreshed with the modified extinction results of 
those layers. 

4. 5.2. 2 Specific Optical Properties Retrieval Procedure for Layers above 20 km 

Layers identified above 20 km are, by GLAS definition, stratospheric aerosol layers or PSCs (which 
are treated as special aerosols). Because the molecular and particulate lidar signal strengths are 
generally low in this region, all profiles in the 20 second block are averaged together before 
processing initiates to sufficiently increase the signal-to-noise ratio. If a layer is found above 20 
km, the top height of the highest 4-second cloud or 4-second aerosol layer below 20 km in the 
current 20-second block is determined. This will then be used to determine the thickness of the 
clear-sky zone below the “above 20 km” layer for possible use in the integrated ratio technique. 

The particulate transmission is initialized to 1 .00 above any first layer found. The molecular 
transmission profile is calculated starting at 41 km altitude where the lidar profile starts. This is, for 
all processing purposes, the top of the atmosphere. The appropriate effective S ratio (either from 
the integrated ratio technique or a default look-up) is calculated for the layer. The multiple 
scattering factor estimate is pulled from the appropriate look-up table, as described in Section 3.7. 
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Detector saturation of the return lidar signal, if any, has to be handled appropriately in the 
processing. Detector saturation would rarely occur in layers found above 20 km, but would occur 
occasionally in thick clouds lower in the troposphere. Each 5-Hertz signal profile bin will have three 
potential states: 1) unsaturated, 2) saturated and no 1064 nm substitution done, and 3) saturated 
with 1064 substitution performed. These three states are determined from two accompanying 
flags: the 5-Hertz saturation flag profile plus the substitution switch. While averaging the lidar 
signal profiles to 4-second and 20-second profiles, a new saturation flag profile is developed. A 
saturation threshold is set (1 0 percent). If less than the threshold percent of the bins being 
averaged are in state 2 or 3, then the saturation flag is set to 0. If only the percent of bins that are 
in state 3 is greater than or equal to the threshold, then the saturation flag is set to 1 . If only the 
percent of bins that are in state 2 is greater than or equal to the threshold, then the saturation flag 
equals 2. If the percent of bins in both states 2 and 3 are greater than or equal to the threshold, 
then the saturation flag is represented by the one with the highest percentage: state 2 is set to 2, 
state 3 is set to 1 . Before processing each layer, perform the following saturation decision tree: 

1 ) Does layer have four or more bins with saturated flag equal to 2? 

> If yes, do not process this or any further layer in current profile 

> If no, does layer have more than 0 bins with saturation flag equal to 2? 

> If yes, use Integrated Ratio Technique subroutine output of T ' 2 at bottom 

of current layer and do not process current layer except to estimate real 
optical depth. All subsequent layers in current profile can be processed 
normally. If Integrated Ratio Technique subroutine could not be used, do 
not process current or any subsequent layer in this profile. 

> If no, proceed 

2) Does layer have three or more bins with saturated flag equal to 1 ? 

> If yes, at the quality flag step at the end of the error processing, bump all 

quality flags of current layer up 2 categories higher 

> If no, does layer have one to two bins with saturated flag equal to 1? 

> If yes, at the quality flag step at the end of the error processing, bump all 
quality flags of current layer up 1 category higher 

> If no, proceed normally (no saturation) 

The particulate transmission calculated at the bottom of the current layer will be used as the initial 
particulate transmission for any lower layer. This transmission is checked to make sure it remains 
inside threshold boundaries before moving to any subsequent layers below. As a final step, the 
optical properties retrievals for any 20-second layer above 20 km are replicated and parceled out to 
the five 4-second standard profiles contained in the 20-second block which are used as the final 
aerosol product resolution. 

4. 5. 2. 3 Specific Optical Properties Retrieval Procedure for Layers below 20 km 

Layers below 20 km are processed much like layers above 20 km. The generic algorithm in 
Section 3.5 is applied and the aerosol S-ratio default assignments are from elevated tropospheric 
or PBL look-up tables. The same detector saturation decision matrix applies that was used above 
20 km. PBL layers are processed down to just above the Earth’s surface. The differences are: 

1 . The starting particulate transmission is the transmission left over from the bottom of the 
last layer above 20 km. 
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2. Data is processed at 4-second resolution, so there are 5 profiles in the 20-second block. 

3. 4-second cloud layers are processed in vertical sequence with the aerosol layers in order 
to accurately calculate the transmission through the clouds. Clouds have their own default 
S ratio equation based on temperature of the cloud. 

4. The PBL can sometimes become a complex layer type of both cloud and aerosol in which 
special rules apply. 

5. 

The PBL is the only layer where cloud and aerosol layers are allowed to overlap. This stems from 
the fact that PBLs are frequently cloud-capped at the top of the inversion. Processing the PBL is 
troublesome in these complex situations. The following rules guide the updated GLAS PBL 
processing: 

A) Do not optically process any PBL which has its cloudy flag turned on (meaning PBL top is 
cloud-covered). 

B) Do not process if a 4-second cloud bottom is at or below the PBL top 

C) Do not process if the current latitude is south of 65S and the PBL quality flag is set to 1 . 
(This eliminates processing of bogus PBLs sometimes found in the Antarctic.) 

A recent modification in the way the GLAS layer search algorithm treats clouds fully imbedded 
inside a PBL (such as cumulus) has allowed for improved aerosol optical processing. In this 
situation, the PBL is reclassified as an elevated aerosol layer whose bottom is the cloud top. 
Aerosol optical processing can then proceed to the embedded cloud top. In earlier versions, the 
PBL above the cloud would not be processed. 

Recent versions of the GLAS GLA1 1 standard product have added the total column aerosol optical 
depth and associated use flag. This optical depth is calculated at night by summing the optical 
depth calculated in the free troposphere with any from aerosol layers processed below the defined 
bottom of the free troposphere. In the sunlit portion of the orbit, the total column aerosol optical 
depth is calculated by only summing the optical depth of all aerosol layers processed. The use flag 
indicates night or day processing plus the layer processing status, such as when the total column 
OD is incomplete because it contains a layer which could not be processed. 

Once the 4-second average aerosol processing has finished, then the 1-sec average cloud 
processing commences for the profile below 20 km using the results of the aerosol analysis to 
determine boundary conditions. 

4.5.2. 4 Specific Column Optical Depth Retrieval Procedure from 1064 nm Surface Return 

Operationally, layer optical depth retrievals from GLAS are limited to the 532 nm atmospheric 
channel. This channel was designed to have the best signal-to-noise and calibration and (through 
a forward lidar inversion algorithm) produce very reasonable (~ 30% error) optical depth analysis of 
every atmospheric particulate layer down to the attenuation of the signal (around 3 optical depth). 
Unfortunately, this channel produced quality profiles for only the Laser 2a (October-November 
2003) period and the first half of the Laser 2b ( February-March 2004) period because of 
deteriorating laser energy for 532 nm in the succeeding Laser 2 and 3 periods. The 532 channel 
was not turned on for the short-lived Laser 1 period. 

Attempts to use the other atmospheric channel at 1064 nm for optical depth retrievals have so far 
been difficult, subject to a noisy signal and an electronic droop effect after strong signals. The 
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1064 nm channel is strong enough for significant layer location detection, but will generally miss 
thin cirrus and aerosol layers (optical depths less than 0.08 - 0.10). 

Fortunately, in addition to the atmospheric scattering, the GLAS measurement includes a precise, 
15 cm resolution, acquisition of the surface waveform at 1064 nm as described in Zwally (2002). In 
fact, the primary science of GLAS involves the use of this waveform for accurate surface altimetry 
work and the fact that both this waveform and the atmospheric profile channels are on the same 
satellite makes GLAS unique. One very useful data product from the waveform is the integrated 
pulse energy from the surface. This received signal by the lidar is a function of the surface 
reflectance and atmospheric transmission. I f one k nows a pr iori or can model the surface 
reflectance with enough precision, then the ratio of the calculated or apparent reflectivity 
(computed from the ratio of the received surface energy to the transmitted laser energy) to the true 
surface reflectivity will be related to the 1064 nm total column atmospheric transmission and thus to 
total column optical depth at 1064 nm. 

In this document we refer to an ocean model of surface reflectance as a function of wind speed 
described and tested by Lancaster (2005) with GLAS data that has shown enough precision to use 
in this approach. A new operational GLAS 1064 nm total optical depth product (on GLA11) is now 
being produced. This new product expands optical depth retrievals beyond the restricted 532 nm 
analysis to cover the Laser 1 period and most of Laser 2 and 3 periods whenever the satellite is 
over ocean and a n on-saturated surface return is detected. In the final release, we have also 
added total column 1064 optical depth over ice sheets by using an assumed, fixed value for the 
surface reflectance of 0.82. While this is indeed the average surface reflectance over ice sheets, it 
does vary considerably between 0.6 and 0.95 and thus the optical depth retrievals over ice sheets 
will be in error by 20-30% at times. 

Details of the algorithm follow. The ocean surface reflectance model we have chosen to use has 
its beginnings with work from Cox and Monk (1954) and Monahan and O'Muircheartaigh (1980). 
At 1064 nm wavelength, ocean reflectance consists predominantly of Fresnel reflection plus a 
small contribution from scattering from whitecaps and sea foam. The ocean surface reflectance 
(R) can be written as: 

(4.5.8) R = (l-W)R s + WR f , 

where R s is the Fresnel reflectance from the surface, Rf is the reflection due to whitecaps, and W is 
the fraction of the surface covered by whitecaps. R s is, in part, a function of the variance of the 
distribution of wave slopes, which is a function of wind speed. The fractional coverage of white 
caps is also a function of wind speed. For details of this equation, refer to Lancaster (2005). The 
updated versions of the GLAS atmospheric data products contain meteorological data interpolated 
from the National Center of Environmental Prediction (NCEP) gridded data set for use in Global 
Climatic Model initialization and contain surface wind speed for every second of orbit track. The 
resultant R from equation 4.5.8 contains no atmospheric attenuation affects from Rayleigh or 
particulate (clouds and aerosol) scattering and is described as the ‘pristine’ surface reflectance one 
would retrieve from a satellite lidar given a known wind speed if there was no atmosphere. Valid 
surface reflectance results are generally limited to values less than 1 .5 because of the instability in 
surface reflectance under calm wind conditions. 
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GLAS Attenuated Backscatter for 10—04—2003 532nm 
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Figure 4.5.4 A Saharan dust event on October 4, 2003 is shown in an image of GLAS backscatter 
profiles in the top plot. The lower plot shows optical depth retrievals from the new 1064 nm surface 
reflectance algorithm (red), 1064 nm GLAS lidar backscatter inversion (blue), and 532 nm GLAS 
standard product lidar backscatter inversion (green). The plot shows the inefficiencies of the 1064 
atmospheric channel to retrieve optical depth and the reasonable correlation between the surface 
reflectance method and the 532 standard optical depth product for dust particulates. 

The GLAS parameter of interest which retrieves the mix of the integrated pulse reflectance from 
the surface and atmospheric attenuation is located in the standard data product GLA05 under the 
name i_reflctUncorr and is calculated at full resolution (40Hz or 175 meters horizontal). Three 
corrections to this parameter must be made before comparing to the pristine surface reflectance for 
particulate optical depth: 

(4.5.9) R G =(R i C b )/(cos(0)T>) 


where Rg is the resultant corrected GLAS reflectance, Ri is the initial ijeflctUncorr, Cb is a boresite 
calibration factor which periodically changes with time, 0 is the tilt angle of the lidar with respect to 
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Backscatter Cross Section (1/m-sr) 


nadir viewing (normally 0.1 but can reach 5.0), and T 2 is the mean molecular two-way 
transmission for the entire atmospheric column at 1064 nm (-0.9853). Both the correction for the 
tilt angle and the molecular transmission are very minor. The relationship between the corrected 
observed GLAS reflectance (Rg) and the modeled pristine ocean reflectance (R) at the observed 
wind is described below: 

(4.5.10) R g =Rq- 2t 

where x is the optical depth of the particulates in the atmospheric column. Solving for x results in 
the equation: 


(4.5.11) r = -hn{R G /R) 

which would be valid for all conditions where the GLAS surface waveform is not saturated, where a 
surface signal is not extinguished by overlying clouds or aerosol, and where the ocean is not 
approaching a windless surface. 


0.0 0.1 0.2 0.3 0.4 0.5 



Figure 4.5.5 Total column aerosol optical depth at 1064 nm based on modeled ocean surface 
reflectivity for the laser 3J operational period. 
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4.5.3 Interpreting the Output 

The output of the optical processing algorithm is separated into two standard product packages, 
GLA1 0 and GLA1 1 . GLA1 0 focuses on the output of vertical profiles: cloud and aerosol 
backscatter cross section and cloud and aerosol extinction cross section. All layer locations are 
referenced from sea level. GLA11 focuses on particulate optical depths (cloud and aerosol). All 
extinction and optical depth values have been corrected for multiple scattering. Polar stratospheric 
clouds in both packages are part of the aerosol category. A list of the key scientific output for 
GLA10 follows: 

14 532 nm aerosol backscatter cross section, 41 to -1 km above mean sea level , 4 second 
averages (zero outside layers) 

15 Aerosol backscatter layer use and quality flags at 4 second resolution, 1 each per layer, 9 
layers, (15 where not processed) [use flag will stipulate saturation state] 

16 532 nm cloud backscatter cross section, 20.4 to -1 km above mean sea level, 1 second 
averages (zero outside layers) 

1 7 Cloud backscatter layer use and quality flags at 1 Hz, 1 each per layer, 1 0 layers, (Value of 
15 where not processed) [use flag will stipulate saturation state] 

18 532 nm aerosol extinction cross section, corrected for multiple scattering, 41 to -1 km 
above mean sea level, 4 second averages (invalid outside layers) 

19 Aerosol extinction layer use and quality flags at 4 second resolution, 1 each per layer, 9 
layers, (15 where not processed) [use flag will stipulate aerosol type] 

20 Calculated aerosol true extinction to backscatter ratios at 4 second resolution, 1 per layer, 

9 layers, (invalid where not processed) [PSC ratio, if layer is PSC] 

21 Default aerosol true extinction to backscatter ratios at 4 second resolution, 1 per layer, 9 
layers, (invalid where not processed) [PSC ratio, if layer is PSC] 

22 Flag stipulating which extinction to backscatter ratio was used in algorithm at 4 second 
resolution, 1 per layer, 9 layers, (invalid where not processed) [PSC ratio, if layer is PSC] 

23 532 nm cloud extinction cross section, corrected for multiple scattering, 20.4 to -1 km 
above mean sea level at 1 Hz, (invalid outside layers) 

24 Cloud extinction layer use and quality flags at 1 Hz, 1 each per layer, 1 0 layers, (1 5 where 
not processed) [use flag will stipulate cloud type] 

25 Calculated cloud true extinction to backscatter ratios at 1 Hz, 1 per layer, 10 layers, (invalid 
where not processed) 

26 Default cloud true extinction to backscatter ratios at 1 Hz, 1 per layer, 10 layers, (invalid 
where not processed) 

27 Flag stipulating which extinction to backscatter ratio was used in algorithm at 1 Hz, 1 per 
layer, 10 layers, (invalid where not processed) 

28 Medium resolution cloud top heights for layers which were selected for optical processing 
at 1 Hz, 1 per layer, 10 layers (invalid where not detected or used) 

29 Medium resolution cloud bottom heights for layers which were selected for optical 
processing at 1 Hz, 1 per layer, 10 layers (invalid where not detected or used) 

30 Medium resolution processed ground detection height at 1 Hz, 1 per profile (invalid where 
not processed) 

31 Low resolution aerosol layer top heights for layers which were selected for optical 
processing at 4 second resolution, 1 per layer, 9 layers, including the planetary boundary 
layer and PSC (invalid where not detected or used) 
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32 Low resolution aerosol layer bottom heights for layers which were selected for optical 
processing at 4 second resolution, 1 per layer, 9 layers, including the planetary boundary 
layer and PSC (invalid where not detected or used) 

33 Low resolution processed ground detection height at 4 second resolution, 1 per profile, 
(invalid where not processed) 

34 Layer top and bottom temperature (C) for all layers sensed 

35 Layer top and bottom pressure (mb) for all layers sensed 

36 Layer top and bottom relative humidity (%) for all layers sensed 

37 Surface meteorological condition (temperature, pressure, relative humidity, and wind) 

38 Bottom height of cloud-free troposphere 

Items 1 through 14 and 25 are calculated by the optical properties algorithm. Items 15 through 24 

are taken from GLA09 and GLA08 particulate boundaries output, but modified to suit the rules 

listed in section 4.5.2 so that only cloud and/or aerosol layers processed optically will show up in 

this data set’s layer locations. 

A list of the key output for GLA1 1 follows: 

1 . 532 nm cloud optical depth, corrected for multiple scattering, at 1 Hz, 1 per layer, 10 layers, 
(invalid where not processed) 

2. Cloud optical depth use and quality flags at 1 Hz, 1 each per layer, 10 layers, (Value of 15 
where not processed) [use flag will stipulate cloud type] 

3. 532 nm elevated aerosol optical depth, corrected for multiple scattering, at 4 second resolution, 
1 per layer, 8 layers, (invalid where not processed) 

4. Elevated aerosol optical depth use and quality flags at 4 second resolution, 1 each per layer, 8 
layers, (15 where not processed) [use flag will stipulate aerosol type, including PSC] 

5. 532 nm planetary boundary layer aerosol optical depth, corrected for multiple scattering, at 4 
second resolution, 1 per layer, 1 layer, (invalid where not processed) 

6. Planetary boundary layer aerosol optical depth use and quality flags at 4 second resolution, 1 
each per layer, 1 layer, (15 where not processed) [use flag will stipulate aerosol type] 

7. Medium resolution 532 nm cloud top heights for layers which were selected for optical 
processing at 1 Hz, 1 per layer, 10 layers, (invalid where not detected or used) 

8. Medium resolution 532 nm cloud bottom heights for layers which were selected for optical 
processing at 1 Hz, 1 per layer, 10 layers, invalid where not detected or used) 

9. Medium resolution 532 nm processed ground detection height at 1 Hz, 1 per profile, (invalid 
where not processed) 

10. Low resolution 532 nm elevated aerosol layer (including PSC) top heights for layers which 
were selected for optical processing at 4 second resolution, 1 per layer, 8 layers (invalid where 
not detected or used) 

1 1 . Low resolution 532 nm elevated aerosol layer (including PSC) bottom heights for layers which 
were selected for optical processing at 4 seconds, 1 per layer, 8 layers (invalid where not 
detected or used) 

12. Low resolution 532 nm processed ground detection height at 4 second resolution, 1 per profile, 
(invalid where not processed) 

13. Low resolution 532 nm planetary boundary layer height at 4 seconds, 1 per profile, (invalid 
where not processed) 
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14. Cloud multiple scattering coefficients used at 1 Hz, 1 per layer, 10 layers, (invalid where not 
processed) 

15. Aerosol multiple scattering coefficients used at 4 second resolution, 1 per layer, 9 layers, 
(invalid where not processed) [including PSC and PBL aerosol] 

16. Multiple scattering effect warning flag at 1 Hz, 1 per profile (Value of 15 where not processed) 

17. Estimated Range (Altimetry) Delay (millimeters) at 1 Hz, 1 per profile taken from last layer 
analyzed 

18. Particle size estimate used to calculate warning flag and range delay at 1 Hz, 1 per profile 
taken from last layer analyzed 

19. Range (Altimetry) Delay Uncertainty (millimeters) at 1 Hz, 1 per profile taken from last layer 
analyzed 

20. Layer top and bottom temperature (C) for all 532 nm layers sensed 

21 . Layer top and bottom pressure (Hp) for all 532 nm layers sensed 

22. Layer top and bottom relative humidity (%) for all 532 nm layers sensed 

23. Medium resolution 1064 nm cloud top heights at 1 Hz, 1 per layer, 10 layers, (invalid where not 
detected or used) 

24. Medium resolution 1064 nm cloud bottom heights at 1 Hz, 1 per layer, 10 layers, invalid where 
not detected or used) 

25. Low resolution 1064 nm aerosol layer top height at 4 second resolution, 1 per layer, 2 layers 
(invalid where not detected or used) 

26. Low resolution 1064 nm aerosol layer bottom height at 4 seconds, 1 per layer, 2 layers (invalid 
where not detected or used) 

27. Layer top and bottom temperature (C) for all 1064 nm layers sensed 

28. Layer top and bottom pressure (mb) for all 1064 nm layers sensed 

29. Layer top and bottom relative humidity (%) for all 1064 nm layers sensed 

30. Surface meteorological condition (temperature, pressure, relative humidity, and wind) 

31. 40 Hz 1064 nm total column optical depth from surface reflectance algorithm 

32. 40 Hz 1064 nm multiple scattering correction factor used 

33. 1 Hz 1064 nm total column optical depth from surface reflectance algorithm 

34. 1 Hz 1064 nm multiple scattering correction factor used 

35. 1 Hz 1064 nm modeled surface reflectance 

36. Total column 532 nm aerosol optical depth at 4 second resolution 

37. Total column 532 nm aerosol optical depth use flag at 4 second resolution 

38. Blowing snow range delay (mm) at 1 Hz 

39. Blowing snow range delay confidence flag 

Items 1 through 6, 14 through 19, and 31 through 39 are calculated by the optical properties 
algorithm. The multiple scattering warning flag and range delay will be based on the height and 
optical depth of the scattering layers and an assumption on the particle size. This is discussed at 
length in section 3.7.2. In general, for clear regions (no cloud or aerosol layers found), the value of 
the multiple scattering warning flag will be zero. The largest value of this flag (14) will occur for 
optically thick layers. Likewise, the range delay will be near zero for clear conditions and a 
maximum for low scattering layers comprised of (assumed) large particles. The details on the 
calculation of the estimated range delay are in section 3.7.2. Items 7 through 1 3 and 20 through 
30 are taken from GLA09 and GLA08 particulate boundaries output. 
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4.5.4 Real Time Error Analysis and Quality Control 

Real time error analysis is performed operationally for every profile analyzed by a simple method 
where an attenuated backscatter profile with error and an S ratio with error are processed through 
the lidar inversion and the optical results with error are compared to the results from the product 
attenuated backscatter profile and the selected product S ratio. 

The backscatter error profiles are produced the same for each horizontal resolution using their 
respective averaged backscatter profiles and standard deviation profiles (see Section 4.5.1. 1 Item 
1). The backscatter profiles are first vertically smoothed using a five bin running average. Next, 
difference profiles are developed by subtracting the smoothed profiles from the averaged 
backscatter profiles, keeping track of sign. Where the difference is zero or positive, the appropriate 
standard deviations are added to the averaged backscatter profile; where it is negative, the 
standard deviations are subtracted from the averaged backscatter profile. The resultant profiles 
are used as the attenuated backscatter error profiles. 

Calculating the S ratio with error is straight forward when conditions allow for the direct calculation 
from the signal loss below the layer. The attenuated backscatter error profile is used instead of the 
regular backscatter profile, resulting in an S ratio with error. Since aerosol S ratios are selected 
from aerosol type cross-reference tables when direct calculation is not feasible, it is harder to 
extract a meaningful error. To calculate the S ratio with error in these situations, the GLAS 
algorithm assigns a five percent deviation which is added to the value pulled from the table. 

Once the attenuated backscatter profiles with error and the S ratios with error are produced, the 
lidar inversion process for the whole atmospheric column is run again with these as inputs, keeping 
track of the new particulate transmission calculations. The resultant extinction, corrected 
backscatter, and optical depths with error are saved temporarily for input into quality flags. 

The GLAS optical processing error analysis culminates into the production of quality flags for 
extinction, corrected backscatter, and optical depth for each layer. The quality flags for all three 
are categories of calculated percent error ranging from the first category (0), which represents 0-5 
% to category 14, which represents 70% and greater. Each category is 5% higher than the last. 
Category 1 5 is saved for the situation where the flag could not be processed. The percent error for 
layer optical depth is simply calculated by subtracting the optical depth with error from the product 
optical depth and dividing by the product optical depth and multiplying by 100. Since the extinction 
and corrected backscatter products are profiles, the percent error is calculated for each vertical bin 
inside the layer, then averaged for the layer. These quality flags are paired with their associated 
use flags per layer. The complete listing of the quality flag and use flag categories can be found in 
the following paragraphs. 

Quality control will be implemented at all stages of the molecular and particulate transmission 
profile development. All input parameters and arrays will be evaluated for quality before being 
used: 

1 . Attenuated backscatter profiles 

• Bad shots will be detected by integration of the lidar signal in the 20 to 40 km height zone. 
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• Lidar bins using 1064 nm backscatter in place of a saturated 532 nm condition will be 
tracked as far as which particulate layer they occur in. 

• Calibration constants which fall outside an expected range will be flagged. 

2. Cloud and aerosol layer detection 

• The layers will be screened so they don’t overlap or become embedded except in the PBL. 

• Visual screening with imagery will occur to make sure layers are labeled ‘cloud’ or ‘aerosol’ 
or ‘polar stratospheric cloud’ correctly. 

3. Molecular backscatter 

• Backscatter calculations from MET data will be monitored to make sure they fall within 
expected boundaries based on atmospheric standards. 

• If MET data are missing or bad, they will be defaulted to atmospheric standards. 

4. Extinction to backscatter ratios and multiple scattering factors 

• The accuracy of these input parameters are at times uncertain, especially for cirrus clouds, 
making this a limitation in the algorithm. 

• Calculations of these parameters in level 2 processing involve a decision matrix look up 
table, which will restrict these parameters to within theoretical and observed limits. If 
atmospheric conditions are favorable, S p will be calculated for thin clouds, elevated 
aerosols, and PSC’s, then compared to matrixed values. If more accurate calculations 
come out of level 3 processing, these will be used to re-process level 2 products. 

As the transmission profiles are processed, the transmission calculations will be tested for out-of- 
bounds situations such as increasing transmission with range or large negative transmission. 
Quality flags will be produced for each particulate layer or profile to help pinpoint how many and 
which output parameters are suspect. This information will be transferred to each of the individual 
output parameter’s quality flags by the following algorithm: 

1 . After calculating optical inversion for layer (eq. 3.5.1 0, etc), re-calculate with a S ratio with 
error (Sp_err) and a signal profile with error (sig_err) as inputs to determine percent error 
for each of the following parameters: 

• Percent backscatter error profile: | Pp-Pp_err| / pp 

• Percent extinction error profile: |ap-ap_err| / ap 

• Percent optical depth error for layer: |xp-xp_err| / xp 

2. Layer quality flags will be set 0-1 5 based on layer averages of the above two percent error 
profiles plus optical depth percent error as follows: 


FLAG 

%ERROR 

FLAG 

%ERROR 

0 

0-5 

12 

60-65 

1 

5-10 

13 

65-70 

2 

10-15 

14 

70 and greater 

3 

15-20 

15 

could not process 

4 

20-25 



5 

25-30 



6 

30-35 



7 

35-40 



8 

40-45 



9 

45-50 
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10 

11 


50-55 

55-60 


3. Layer use flags are designated as follows: 

a) for backscatter cross section, the use flag gives saturation status as follows: 

FLAG SATURATION STATUS 

0 no saturation detected 

1 one or two bins were saturated with 1 064 nm conversion performed 

2 at least three bins were saturated with 1064 nm conversion performed 

3 at least one but less than four bins were saturated with no conversion 
performed 

4 four or more bins were saturated with no conversion performed 

15 invalid 

b) for extinction cross section and layer optical depth, the use flag designates layer type 
category as follows: 

Stratospheric Aerosol: {Layers 1-3} 

FLAG CATEGORY 

00 Generic default 

12 STRATO aerosol (any non-PSC layer whose top is > tropopause 

1 3 PSC type I (PSC with rh less than or equal to 95%) 

14 PSC type II (PSC with rh greater than 95%) 

15 invalid 

Tropospheric Aerosol: {Layers 4-9} 

FLAG CATEGORY 

00 Generic default 

01 Sulfate+carbon 

02 Carbon 

03 Salt+dust 

04 Salt 

05 Sulfate 

06 Dust+carbon 

07 Salt+dust+sulfate 

08 Salt+carbon 

09 Salt+sulfate 

10 Dust 

1 1 Salt+dust+carbon 

12 Dust+sulfate 

13 Salt+carbon+sulfate 

14 Dust+carbon+sulfate 

15 Mix of all types 

Cloud: {based on average cloud temperature, water cloud is warmer than -1 3 C} 

FLAG CATEGORY 

00 less than or equal to -75.0 C 

01 -75.0 through -68.5 
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02 -68.5 through -62.0 

03 -62.0 through -55.5 

04 -55.5 through -49.0 

05 -49.0 through -32.5 

06 -32.5 through -26.0 

07 -26.0 through -19.5 

08 -19.5 through -13.0 

09 -13.0 through -6.5 

10 -6.5 through 0.0 

11 0.0 through 6.5 

12 6.5 through 13.0 

13 13.0 through 19.5 

14 greater than 19.5 C 

15 invalid 
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5 Mitigating Multiple Scattering Induced Ranging Errors 

It has been calculated that the effects of multiple scattering from cloud and 
aerosol will introduce significant errors for precision surface altimetry. These 
results are presented in detail by Duda et al. (1999, a and b and available from 
the GLAS ftp site). The pulse spreading from multiple scattering will tend to 
introduce a positive bias to the range determination. The magnitude of the 
effect can be considerable under certain atmospheric conditions, ranging to 
larger than 1 meter for a single pulse depending on conditions. Since cloud 
cover varies seasonally and year to year, Duda et al. show that if uncorrected, 
the multiple scattering effect would introduce significant errors for the GLAS 
surface altimetry yearly analyses. The magnitude of the range delay effect is a 
function of the scattering layer height, optical depth and the physical 
characteristics of the scattering particles (size, shape, composition). What was 
found from analysis of GLAS data is that over Antarctica, the major cause of 
range delay was from blowing snow layers only 100 - 200 m or so thick. These 
layers are in contact with the ground and can have optical depths approaching 
1.0. Figure 5.1 below shows an example of the dramatic ranging errors (nearing 
10 m) that can occur due to multiple scattering of photons through blowing 
snow layers. 

Application of the atmospheric channel of GLAS to perform an analytic 
correction to the multiple scattering induced ranging error has been developed. 
In section 3.7.2 we present a detailed discussion of the approach which is 
based on the creation of a ranging error table for many different atmospheric 
conditions. The determining factors are the cloud height range and optical 
thickness plus an assumption of cloud particle size. The factors are essentially 
the same as those to be used for the generation of the correction factors for the 
influence of multiple scattering on cloud and aerosol cross sections and optical 
thickness as also described in section 3.7. As described in Duda et al., an 
estimate of the magnitude of the pulse spreading error on the surface is 
computed based on a centroid analysis of a flat, normal surface. This 
information can then be used by the altimetry processing to eliminate shots that 
are likely to be severely affected by multiple scattering. GLA1 1 will have 3 items 
in its output that are related to the multiple scattering induced range delay. They 
are 1) Multiple scattering warning flag, 2) Particle size used in multiple 
scattering computation and 3) Calculated range-to-surface delay. They are 
described in detail in section 3.7.2 and listed in section 4.5.3. 
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Figure 5.1 Blowing snow layer detected by ICESat over southern Greenland on 23 October 2003. 
(b) shows the 532 nm attenuated backscatter along the segment shown on the map (red line on 
map) and (a) shows the difference in ICESat measured elevations between 23 October 2003 and 
1 1 October 2004 along the same track (white '+’ signs). Also shown in (a) is the optical depth of the 
blowing snow layer (red line and axis on right). 
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6 Future Research 


Unfortunately, the instrument problems encountered by GLAS have severely 
curtailed the number of things that can be done with the atmospheric data. 
Nevertheless, the level II data products discussed in this document form the 
basis for quite a few scientific studies over the last 5-6 years. To name a few, 
the data have been used to study polar clouds (Palm et al., 2010), Boundary 
layer height and verification of GCM models (Palm et al., 2005), blowing snow 
over Antarctica (Palm et al., 2011), tropospheric cirrus (Dessler et al., 2006), and 
Polar Stratospheric Clouds (Palm et al., 2005; Mueller et al., 2008). Future 
research activities could include a further exploitation of the global cloud data 
set, and the vertical distribution of aerosol. Another area that could be mined 
from the data is an extension of the PBL height analysis to include a derivation 
of lifting condensation level (LCL), derived from the output of GLA08. The more 
difficult parameters to obtain accurately from the lidar data are the optical depth 
and extinction cross sections for aerosol and cloud. It is expected that the 
accuracy and applicability of these can be significantly increased through Level 
III products and post processing. The two areas requiring further work for this 
are the use of data other than the lidar profile signal and improvements in 
multiple scattering corrections. 

For the cloud analysis, a desirable input would be simultaneous IR radiance 
measurements. With IR radiance obtained in sufficiently close time with the lidar 
profile it is possible to solve for the vertical profiles of IR absorption cross 
section (Spinhirne et al., 1990). Simultaneous IR radiance values are available 
for a large fraction of the GLAS observations. During the mission lifetime, there 
were over 20 spectral imagers with thermal IR channels in orbit. Since GLAS 
has a precessing orbit, the GLAS measurements will be within the swath width 
of the MODIS imagers for about two months of the year for example. The 
combination of the GLAS data with IR data could be a fruitful research topic for 
future research. An additional improvement of the cloud retrieval from GLAS 
data alone may also be possible from research and modeling on using the 
molecular and surface signals under thin cloud layers to improve optical depth 
calculations. The most significant improvement for cloud retrieval will likely 
come from research on the best approach for the multiple scattering correction. 
One area to explore is the use of the below ground multiple scattering tail that 
has been observed by the GLAS ranging channel for a direct measure of the 
multiple scatter factor leading to improvement of correction tables. Elevation 
retrievals through blowing snow layers should be compared with the repeat 
track elevation values for clear conditions to compare with model calculations 
and the range delays calculated and reported on the atmospheric products 
which could lead to the development of better methods to correct for the 
multiple scattering induced errors for surface ranging. 
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For aerosol optical thickness and extinction cross section, multiple scattering 
corrections are less of an issue. The largest uncertainties would result likely 
from the value of extinction to backscatter ratio that is applied for the retrievals. 
An important factor for improving the retrievals for day time observations is to 
make use of the 532 and 1064 nm solar background signals. From these data 
alone, over oceans an optical thickness for aerosol could be obtained directly in 
the manner that AVHRR data are now used. Future research is needed to model 
the best approach for incorporating the solar background signals with the lidar 
return profiles. Yang et al., 2008 have already demonstrated the use of GLAS 
532 background data to retrieve thick cloud optical depth. In addition, the GLAS 
aerosol profiles can be combined with many other sensor data and retrievals. 
One example would be with AVHRR and MODIS aerosol retrievals. Again the 
precessing orbit of GLAS will provide large amount of coincident data that can 
be used to improve extinction to backscatter look up tables for nighttime and 
other non-coincident GLAS observations. An especially important combination 
will be GLAS aerosol profiles with TOMS aerosol retrievals. Currently TOMS 
data are applied to retrieve absorbing aerosol in the atmosphere, but an 
assumption on the height profile of the aerosol is needed. For the large amount 
of coincident data with TOMS expected from the full GLAS mission, future 
research will enable improvements in the TOMS and GLAS aerosol data results. 


128 



7 References 


Ackermann, J., 1998: The Extinction-to-Backscatter Ratio of Tropospheric 
Aerosol: A Numerical Study, J. Atmos. Oceanic Technol., 15, 1043-1050. 

Ansmann, A., U. Wandinger, M. Riebesell, C. Weitkamp and W. Michaelis, 1992 
Independent measurements of extinction and backscatter profiles in cirrus 
clouds by using a combined Raman elastic-backscatter lidar, Appl. Opt., 31, 
7113-7131. 

Ansmann, A., J. Bosenberg, G. Brogniez, S. Elouragini, P. Flamant, K. Klapheck, 
H. Linn, L. Menenger, W. Michaelis, M. Riebesell, C. Senff, P. Thro, U. 
Wandinger and C. Weitkamp, 1993: Lidar Network Observations of Cirrus 
Morphological and Scattering Properties during the International Cirrus 
Experiment 1989: The 18 October 1989 Case Study and Statistical Analysis, J. 
Appl. Meteor., 32, 1608-1622. 

Boers, R, E. W. Eloranta and R. L. Coulter, 1984: Lidar observations of mixed 
layer dynamics: tests of parameterized entrainment models of mixed layer 
growth rate. J. Clim. Appl. Meteor, 23, 247-266. 

Boers, R. and E. W. Eloranta, 1986: Lidar measurements of the atmospheric 
entrainment zone and the potential temperature jump across the top of the 
mixed layer. Bound. Layer Meteor., 34, 357-375. 

Boers, R., S.H. Melfi and S.P. Palm, 1991: Cold Air outbreak during GALE: Lidar 
observations and modeling of boundary layer dynamics. Mon. Wea. Rev., 119, 
1132-1150. 

Congeduti, F., J. DeLuisi, T. DeFoor, and L. Thomason, 1998: Optical extinction 
properties of volcanic stratospheric aerosol derived from ground-based lidar 
and Sun photometer measurements, J. Geophys. Res., 103, 13893-13902. 

Cox, C., and W. Munk, 1954: Measurement of the roughness of the sea surface 
from photo-graphs of the Sun’s glitter, J. Opt. Soc. Am., 44, 838-850. 

Dessler, A.E., S.P. Palm, W.D. Hart, and J. D. Spinhirne, Tropopause-level thin 
cirrus coverage revealed by ICESat/Geoscience Laser Altimeter System, J. 
Geophys. Res., Ill, 2006. 

Duda, D. P., J. D. Spinhirne and E. W. Eloranta, 1999: Atmospheric Multiple 
Scattering Effects on Altimetry. Part I: Calculation of Single Pulse Bias. 
Submitted February, 1999 

Duda, D. P., J. D. Spinhirne and E. W. Eloranta, 1999: Atmospheric Multiple 
Scattering Effects on Altimetry. Part II: Cloud Climatology Analysis of Expected 
Seasonal and Interannual Surface Altitude Errors. Submitted February, 1999 


129 



Eloranta, E. W., R. E. Kuehn and R. E. Holz, “Cirrus Cloud Backscatter Phase 
Functions Measured with the University of Wisconsin High Spectral Resolution 
Lidar” 10th Conference on Atmospheric Radiation, preprint, AMS, Madison, 
Wisconsin, 28 June-2 July, 1999. 

Elouragini, S., 1995: Useful algorithms to derive the optical properties of clouds 
from a backscatter lidar return, J. Modern. Opt., 42, 1439-1446 

Elouragini, S. and P. H. Flamant, 1996: Iterative method to determine an 
averaged backscatter-to-extinction ratio in cirrus clouds, Appl. Opt., 35, 1512- 
1518. 

Gobbi, G. P., 1995: Lidar estimation of stratospheric aerosol properties: Surface, 
volume, and extinction to backscatter ratio, J. Geophys. Res., 100, 11219- 
11235. 

Grund, C. J. and E. W. Eloranta, 1990: The 27-28 October 1986 FIRE IFO Cirrus 
Case Study: Cloud Optical Properties Determined by High Spectral Resolution 
Lidar, Mon. Wea. Rev., 118, 2344-2355. 

Grund, C.J. and E. W. Eloranta, 1991: The University of Wisconsin High Spectral 
Resolution Lidar, Optical Engineering, 30, 6-12. 

Hlavka, D. L., J. D. Spinhirne, and J. R. Campbell, 1998: “Aerosol Analysis 
Techniques and Results form Micro Pulse Lidar”, 19th International Laser Radar 
Conference, Annapolis, MD, July 6 - 10, 1998. 

Holben BN, Eck TF, Slutsker I, Tanre D, Buis JP, Setzer A, Vermote E, Reagan 
JA, Kaufman YJ, NakajimaT, Lavenu F, Jankowiak I, Smirnov A, 1998: 
AERONET - A federated instrument network and data archive for aerosol 
characterization. Remote Sensing of Environment 66: (1) 1-16. 

IPCC, 1995: Climate Change 1994. Cambridge U. Press. 

Iqbal, M., An Introduction to Solar Radiation, Academic Press, New York, NY, 
1983. 

Lancaster, R. S., J. D. Spinhirne, and S. P. Palm, 2005: Laser pulse reflectance 
of the ocean surface from the GLAS satellite lidar, Geophys. Res. Lett., 32, 
L22S10, doi:10.1029/2005GL023732. 

Liu, Z., P. Voelger, and N. Sugimoto, 2000: Simulations of the observation of 
clouds and aerosols with the Experimental Lidar in Space Equipment system. 
Appl. Opt., 39, 3120-3137. 

Macke, A., 1993: Scattering of light by polyhedral ice crystals, Appl. Opt., 32, 
2780-2788. 


130 



Macke, A., P. N. Francis, G. M. McFarquhar and S. Kinne, 1998: The role of 
particle shapes and size distributions in the single scattering properties of cirrus 
clouds, J. of Atmos. Sci., 55, 2874-2883. 

Marenco, F., V. Santacesaria, A. F. Bais, D. Balis, A. di Sarra, A. Papayannis, 
and C. Zerefas, 1997: Optical properties of tropospheric aerosols determined by 
lidar and spectrophotometric measurements (Photochemical Activity and Solar 
Ultraviolet Radiation campaign), Appl. Opt., 36, 6875-6886. 

McClatchey, R.A., R.W. Fenn, J.E.A. Selby, F.E. Volz, and J.S. Garing, Optical 
Properties of the Atmosphere (Third Edition), AFCRL-72-0497, 1972. 

McCormick, M. P., H. M. Steele, W Chu, and T. Swissler, 1982: Polar 
Stratospheric Cloud Sightings by SAM II. J. Atmos. Sci., 39, 1387-1397. 

McCormick, M. P., P. Hamill and U. Farrukh, 1985: Characteristics of Polar 
Stratospheric Clouds as observed by SAM II, SAGE, and Lidar, J. Meteor. Soc. 
Japan, 63, 267-276. 

McFarquhar, G. M. and A. J. Heymsfield, 1997: Parameterization of Tropical 
Cirrus Ice Crystal Size Distributions and Implications for Radiative Transfer: 
Results from CEPEX, J. Atmos. Sci., 54, 2187-2200. 

Melfi, S.H., J.D. Spinhirne, S.H. Chou and S. Palm, 1985: Lidar observations of 
vertically organized convection in the planetary boundary layer over the ocean. J. 
Clim. Appl. Meteor., 24, 806-821. 

Mishchenko, M. I., Wielaard, and B. E. Carlson, 1997: T-matrix computations of 
zenith-enhanced lidar backscatter from horizontally oriented ice plates, Geophys. 
Res. Lett., 24, 771-774. 

Monahan, E. C., and I. O’Muircheartaigh, 1980: Optimal power-law description 
of oceanic whitecap coverage dependence on wind speed, J. Phys. Oceanogr., 
10,2094-2099. 

Mueller, K.J., L. Di Girolamo. M. Fromm and S.P. Palm, Stereo observations of 
polar stratospheric clouds. Geophys. Res. Lett., 35, 2008. 

Nicolas, F., Bissonnette, L. R., and P. H. Flamant, 1997: Lidar effective multiple- 
scattering coefficients in cirrus clouds, Appl. Opt., 36, 3458-3468. 

Palm, S.P. and J.D. Spinhirne, 1987: Optimization of boundary layer height 
retrieval. Laser and Optical Remote Sensing of the Atmosphere, Volume 18, 63- 
66, and presented at conference. 

Palm, S.P. and J. Spinhirne, 1998: The detection of Clouds, Aerosol and Marine 
Atmospheric Boundary Layer Characteristics from Simulated GLAS Data. The 
19th International Laser Radar Conference, Annapolis, Md, July 6-10, 1998. 


131 



Palm, S. P., D. Hagan, G. Schwemmer and S.H. Melfi, 1998: Inference of Marine 
Atmospheric Boundary Layer Moisture and Temperature Structure using 
Airborne Lidar and Infrared Radiometer Data, J. Appl. Meteor., 37, 308-324. 

Palm, S. P., A. Benedetti, and J. Spinhirne, Validation of ECMWF global forecast 
model parameters using the Geoscience Laser Altimeter System (GLAS) 
atmospheric channel measurements. Geophys. Res. Lett., 32, 2005. 

Palm, S.P. M. Fromm and J. Spinhirne, Observations of Antarctic Polar 
Stratospheric Clouds by the Geoscience Laser Altimeter System (GLAS), 
Geophys. Res. Lett., 32, 2005. 

Palm, S.P., S. Strey, J.D. Spinhirne, 2010, Influence of arctic sea ice extent on 
polar cloud fraction and vertical structure and implications for regional climate, J. 
Geophys. Res., 115, D21209, doi:10.1029/2010JD013900 

Palm, S.P., Y. Yang, J.D. Spinhirne and A. Marshak, 201 1 , Satellite Remote 
Sensing of Blowing Snow Properties over Antarctica, J. Geophys. Res., 
accepted May 2011. 

Pinnick, R. G., J. M. Rosen and D. J. Hofmann, 1976: Stratospheric Aerosol 
Measurements III: Optical Model Calculations, J. Atmos. Sci., 33, 304-314. 

Pinnick, R. G., S. G. Jennings, P. Chylek, C. Ham and W. T. Grandy, 1983: 
Backscatter and Extinction in Water Clouds. J. Geophys. Res, 88, 6787-6796. 

Platt, C. M. R., 1979: Remote Sounding of High Clouds: I. Calculation of Visible 
and Infrared Optical Properties from Lidar and Radiometer Measurements, J. 
Appl. Meteor., 18, 1130-1143. 

Platt, C. M. R., Reynolds, D. W., Abshire, N. L., 1980: Satellite and lidar 
observations of the albedo, emittance, and optical depth of cirrus compared to 
model calculations. Mon. Wea. Rev., 108, 195-204 

Platt, C. M. R., 1981 : Remote sounding of high clouds. Part III: Monte Carlo 
calculations of multiple-scattered lidar returns, J. Atmos. Sci., 38, 156-167. 

Powell, D. M., J. A. Reagan, M. A. Rubio, W. H. Erxleben, and J. D. Spinhirne, 
2000: ACE-2 Multiple Angle Micro-pulse Lidar Observations from Las Galletas, 
Tenerife, Canary Islands. Tellus B, 52, 651-660. 

Russell, P. B., T. J. Swissler, M. P. McCormick, W. P. Chu, J. M. Livingston and 
T. J. Pepin, 1981 : Satellite and Correlative Measurements of the Stratospheric 
Aerosol. I: An Optical Model for Data Conversions, J. Atmos. Sci., 38, 1279- 
1294. 

Spinhirne, J. D., J. A. Reagan, and B. M. Herman, 1980: Vertical Distribution of 
Aerosol Extinction Cross Section and Inference of Aerosol Imaginary Index in the 
Troposphere by Lidar Technique, J. Appl. Meteor., 19, 426-438. 


132 



Spinhirne, J. D., 1982: Lidar clear atmosphere multiple scattering dependence on receiver 
range, Appl. Opt., 21, 2467-2468. 


Spinhirne, J. D., R. Boers and W. D. Hart, 1989: Cloud Top Liquid Water from 
lidar Observations of Marine Stratocumulus. J. Appl. Meteor., 28, 81-90. 

Spinhirne, J. D. and W. D. Hart, 1990: Cirrus structure and radiative parameters 
from airborne lidar and spectral radiometer observations: the 28 October 1986 
FIRE study. Mon. Wea. Rev, 118, 2329-2343 

Spinhirne, J. D., W. D. Hart, D. L. Hlavka, 1996: Cirrus infrared parameters and 
shortwave reflectance relations from observations, J. of Atmos. Sci., 53, 1438- 
1458. 

Spinhirne, J. D. and S. P. Palm, 1996: Space Based Atmospheric Measurements 
by GLAS. Advances in Atmospheric Remote Sensing with Lidar. Selected 
Papers of the 18th International Laser Radar Conference, Berlin Germany. 213- 
216. 

Spinhirne, J. D., Campbell, J. R., Hlavka, D. L., Ferrare, R. A., Turner, D. D., and 
Flynn, C. J., "Aerosol Retrieval Comparison During the SGP Summer ’98 IOP 
from Multiple Lidar Probing", Poster Abstract, Atmospheric Radiation 
Measurement (ARM) Science Team Meeting, San Antonio, Texas, March 22-26, 
1999. 

Starkov, A. V., and C. Flesia, 1998: Correction of spaceborne lidar signal for 
multiple scattering from high clouds, Proc. of the 19th International Laser Radar 
Conference, July 6-10, Annapolis, MD. 

Vigroux, E., Contribution a I'etude experimentale de I'absorption de I'ozone, Ann. 
Phys., 8, 709-761, 1953. 

Welton, E. J., K. J. Voss, H. R. Gordon, H. Maring, A. Smirnov, B. Holben, B. 
Schmid, J. M. Livingston, P. B. Russell, P. A. Durkee, P. Formenti, and M. O. 
Andreae, 2000: Ground-based Lidar Measurements of Aerosols During ACE-2: 
Instrument Description, Results, and Comparisons with other Ground-based 
and Airborne Measurements. Tellus B, 52, 635-650. 

Wielicki, B. A., B. Barkstrom, E. Harrison, R. Lee, G. Smith, and J. Cooper, 1996: 
Clouds and the Earth’s radiant energy system (CERES): An Earth observing 
system experiment. Bull. Amer. Meteor. Soc., 77 

Wiscombe, W., 1977: Mie scattering calculations: Advances in technique and 
fast, vector-speed computer codes. NCAR Tech. Note TN140+STR. (Edited and 
revised 1996, available at 

ftp://climate.gsfc.nasa.gov/pub/wiscombe/Single_Scatt/Mie_Code/NCARMieRe 

port.pdf 


133 



Yang, Y, A. Marshak, J. Chiu, W. Wiscombe, S.P. Palm, A. Davis, D. 
Spangenberg, L. Nguyen, J. Spinhirne and P. Minnis, Retrievals of thick cloud 
optical depth from the Geoscience Laser Altimeter System (GLAS) by calibration 
of solar background signal, J Atmos Sci, 65, 3513->3527, 2008. 

Zwally, H. J., et al., 2002: Ice, cloud and land elevation satellite’s laser 
measurements of polar ice, atmosphere, ocean and land, J. Geodyn., 34, 405- 
445. 


134 



8 Acronyms 


ACE 

AEROCE 

AERONET 

AIRS 

ARM 

ARMCAS 

ASTEX 

AVHRR 

AVRIS 

CALS 

CAR 

CEPEX 

CRYSTAL 

DAAC 

DEM 

EAL 

EOS 

EOSDIS 

FIRE 

GLAS 

GLOBE 

HSRL 

ISCCP 

LITE 

MAS 


Arctic Clouds Experiment 

Aerosol/Ocean Chemistry Experiment 

Aerosol Robotic Network 

Atmospheric Infrared Sounder 

Atmospheric Radiation Measurement Program 

Arctic Radiation Measurements in Column Atmosphere-surface 
System (beaufort Sea, Alaska, June 1995) 

Atlantic Stratocumulus Transition Experiment (Azores, June 1992) 

Advanced Very High Resolution Radiometer 

Airborne Visible / Infrared Imaging Spectrometer 

Cloud and Aerosol Lidar System 

Cloud Absorption Radiometer 

Central Equatorial Pacific Experiment 

Distributed Active Archive Center 

Digital Elevation Model 

Elevated Aerosol Layer 

Earth Observing System 

EOS Data and Information System 

First ISCCP Regional Experiment 

Geoscience Laser Altimeter System 

Global Backscatter Experiment 

High Spectral Resolution Lidar 

International Satellite Cloud Climatology Project 

Lidar In-space Technology Experiment 

MODIS Airborne Simulator 
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MAST 

MODIS 

NSIDC 

PBL 

SUCCESS 

WINCE 


Monterey Area Ship Tracks Experiment (Monterey California, June 
1994) 

Moderate Resolution Imaging Spectroradiometer 
National Ice and Snow Data Center 
Planetary Boundary Layer 

Subsonic Aircraft Contrail and Cloud effects Special Study (April - 
May, 1996) 

Winter Cloud Experiment 
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